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Abstract 


Readmission following discharge from an initial hospitalization is a key marker 
of quality of health care in the United States. For the most part, readmission 
has been studied among patients with ‘acute’ health conditions, such as pneumo¬ 
nia and heart failu re, with anal yses based on a logistic-Normal generalized linear 


mixed model (jNormand et all 19971 1. Naive application of this model to the study 


of readmission among patients with ‘advanced’ health conditions such as pancre¬ 
atic cancer, however, is problematic because it ignores death as a competing risk. 
A more appropriate analysis is to imbed such a study within the semi-conrpeting 
risks framework. To our knowledge, however, no comprehensive statistical methods 
have been developed for cluster-correlated semi-conrpeting risks data. To resolve this 
gap in the literature we propose a novel hierarchical modeling framework for the 
analysis of cluster-correlated semi-competing risks data that permits parametric or 
non-parametric specifications for a range of components giving analysts substantial 
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flexibility as they consider their own analyses. Estimation and inference is performed 
within the Bayesian paradigm since it facilitates the straightforward characterization 
of (posterior) uncertainty for all model parameters, including hospital-specific ran¬ 
dom effects. Model comparison and choice is performed via the deviance information 
criterion and the log-pseudo marginal likelihood statistic, both of which are based 
on a partially marginalized likelihood. An efficient computational scheme, based on 
the Metropolis-Hastings-Green algorithm, is developed and had been implemented 
in the SemiCompRisks R package. A comprehensive simulation study shows that the 
proposed framework performs very well in a range of data scenarios, and outperforms 
competitor analysis strategies. The proposed framework is motivated by and illus¬ 
trated with an on-going study of the risk of readmission among Medicare beneficiaries 
diagnosed with pancreatic cancer. Using data on n=5,298 patients at J=112 hospi¬ 
tals in the six New England states between 2000-2009, key scientific questions we 
consider include the role of patient-level risk factors on the risk of readmission and 
the extent of variation in risk across hospitals not explained by differences in patient 
case-mix. 

Keywords: Bayesian survival analysis; cluster-correlated data; illness-death models; re¬ 
versible jump Markov chain Monte Carlo; shared frailty; semi-competing risks 
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1 Introduction 


Cancer of the pancreas is one of the most deadly. In 2013, an estimated 38,460 in¬ 
dividuals died from pancreatic c ancer in the Un ited State s making it the fourth most 


prevalent cause of cancer death (jAmerican Cancer Socie ty!, 


20131 ). Unfortunately, since 


there are no effective screening tests for pancreatic cancer, most patients are diagnosed 
at a late stage of the disease, specifically once it has metastasized to other parts of the 
body. As a result, survival is poor w ith 1-year and 5-year mortality rates are of 75% 
and 94%, respectively ([Shin and Cantol. 2012)- In practice, since prognosis is poor and 


mortality rates high, the treatment and management of patients diagnosed with pancre¬ 
atic cancer gener ally f ocuses on palliative care aimed at enhancing quality of end-of-life 
care flPLoS Medicine Editors! . 12012). Such care is expensive, however, with patients diag¬ 
nosed with pan creati c cancer ac cruing an estimated $165,000 in health care costs in their 


last year of life flMar i otto et ah 


201l|). 


Despite the huge costs, there are currently no comprehensive national efforts to monitor 
quality of end-of-life care for pancreatic cancer nor for any of a broad range of other 
‘advanced’ health conditions for which the management of disease focuses on palliative 
care. Outside the context of these conditions, however, there is substantial interest in 
understanding variation in quality of health care. The recent literature, in particular, has 
focused on readmission as a key marker of quality of care, in part becau se it is expensive bu 


also 


Decause it is t 


2011 


rough 


Brooks et al 


2014; 


o f as an ofte n-pre venta blc event flVest et al. 


2010 


Warren et al. 


Stitzenberg et al 


201 5! ). In addition, as the nation’s largest payer 


of health care costs in the United States, the Centers for Medicare and Medicaid Services 
(CMS) uses hospital-specific readmission rates as a central component in two programs: the 
Hospital Inpatient Quality Reporting Program, which requires hospitals to annually report, 
among other measures, readmission rates for pneumonia, heart failure and my ocard ial 
infarction in order to receive a full update to their reimbursement payments fjCMSl . 2013al ): 
and, the Readmission Reduction Program, w hich requires CMS to reduce payments to 
hospitals with excess readmissions fjCMSl . 2013bi). 

Across all of these efforts, investigations of readmission in the literature have invariably 
used a logistic-Normal generalized linear mixed model (LN-GLMM) to analyze patients 
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clustered within hospitals ( Normand et all 1997; 


Ash et al. 


2012 ). While reasonable for 


health conditions with effective treatment options and low mortality, direct application of 
this model to investigate variation in risk of readmission following a diagnosis of pancreatic 
cancer is inappropriate because of the strong force of mortality. Consider, for example, 
n=5,298 Medicare beneficiaries diagnosed with pancreatic cancer at J=112 hospitals in 
six New England states between 2000-2009 and suppose interest lies in understanding 
determinants of readmission 90 days post-discharge. While additional detail is given below, 
we note at the outset that 1,257 patients (24%) died within 30 days of discharge without 
experiencing a readmission event; furthermore, 1,912 patients (36%) died within 90 days 
of discharge without experiencing a readmission event. Naive application of a standard 
LN-GLMM to these data ignores the fact that a substantial portion of the patients are not 
at risk to experience the event of ‘readmission by 90 days’ for much of the timeframe. Such 
an analysis may lead to bias and, if incorporated into existing CMS programs, could have 
a major impact on how hospitals are penalized for poor quality of care. 

In the statistics literature, data that arise from studies in which primary scientific in¬ 
terest lies with some non-terminal event (e.g. readmission) whose observation is subjec 


to & terminal event (e.g. death) are referred to as semi-competing risks data (IFine et al. 


20 011). Broadly, published methods for the analysis of semi-competing risks data can be 


classified into three groups: metho ds that spe c ify dependence between the non-terminal 


and terminal events via a copula ( Fine et al. 


2001 


Peng and Find . 2007; 


Hsieh et al. 


specific frailty ( 

Kncib and Hennerfe 

nd 

, 2008; 

Xu et ah. 2010: Zeng et ah. 

2012 

Han et ah. 

2014; 

Zhang: et al.. 

2014 


Lee et al. 

• 2 

015 

): and, methods based on principal stratifica- 

tion ( 

Zhang and Rubin 

5 

2003; 

Egleston et al. 

2007; Tchetgen Tchetgen. 

2014 

. Common 


to all of these methods, however, is that their development has focused exclusively on 
settings where individual study units are independent. As such, the methods are not de¬ 
sign t o address sci e ntific questions that arise natu rally in the context of cluster-correlated 


data (Diggle et al. 


2002 ; 


Fitzmauricc et al. 


20121). In the context of readmission following 


a diagnosis of pancreatic cancer, such questions include: (i) the investigation of between- 
and within-hospital risks factors for readmission while acknowledging death as a competing 
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force, (ii) characterizing and quantifying between-hospital variation in risk of the terminal 
event not explained by differences in patient case-mix, and (iii) estimating, and quantifying 
uncertainty for, hospital-specific effects, as well as ranking. Furthermore, it is well-known 
that if one is to perform valid inference all potential sources of correlation must be ac¬ 
counted for in the analysis. 

To our knowledge, while the literature on the relat ed co mpeting risks pro 


sider e d methods for cl u ster-correlated da t a set t ings (IKatsahian et al. 


2008 


Gorhne and Hs 


2 


2011 


Zhou et al. 


2012 


Gorfine et al. 


cs prob 

lem has con- 

, 2006; 

Chen et al.. 


2014]), only one paper on 


the anal ysis of clu ster-c orrelated semi-competing risks data has been published. Specifi¬ 
cally, Liquet et alj (2012) recently proposed a multi-state model that incorporated a hospital- 
specific random effect to account for cluster-correlation. Estimation and inference was per¬ 


formed within the frequentist paradigm, based on an integrated likcli 
over the random effect, implemented in the f railtypack R package I 


rood t hat marginalizes 


Rondeau et al. 


20 m 


For our purposes, however, their appro ach is limited in a number of important ways. First, 
the analyses presented in Liq uet et al. (12012 4 permit either a patient-specific frailty to ac¬ 
count for dependence between T\ and T 2 or a hospital-specific random effect to account 
for cluster-correlation but not both simultaneously. Second, the proposed specification as¬ 
sumed that the hospital-specific random effect for the non-terminal event is independent of 
the hospital-specific random effect for the terminal event, precluding a potentially impor¬ 
tant form of dependence. Third, towards understanding variation in risk of readmission, 
the hospital-specific random effects are themselves key parameters of scientific interest and 
not nuisance parameters to be marginalized over. Finally, evaluation of the integrated like¬ 
lihood requires the specification of a parametric distribution for the hospital-specific ran¬ 
dom effects. While estimation and inference for regression parameters is generally robust 
to misspecification of random effects distributions in GLMMs, misspecification is known 


to adv ersely impact the shape of the estimated distribut i on o 


selves (McCulloch et al. 


2011 


McCulloch and Neuhaus, 


2011 


he random effects them- 


Neuhaus and McCulloch, 


2 0111 ). This is particularly important in quality of health care studies where identifying 


a hospital as being in the tail of the distribution can have a substantial impact on their 
evaluation. 
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Towards overcoming these limitations, we develop a novel, comprehensive hierarchical 
multi-state modeling framework for cluster-correlated semi-competing risks data. A key 
feature of the framework, and its implementation, is that it permits either parametric or 
non-parametric specifications for a range of model components, including baseline hazard 
functions and distributions for hospital-specific random effects. This gives analysts substan¬ 
tial flexibility as they consider their own analyses. Estimation and inference is performed 
within the Bayesian paradigm which facilities the straightforward quantification of uncer¬ 
tainty for all model parameters, including hospital-specific random effects and variance 
components. The remainder of this paper paper is organized as follows. Section [2] intro¬ 
duces an on-going study of readmission among patients diagnosed with pancreatic cancer, 
and provides a description of the available Medicare data. Section [3] describes the proposed 
framework, including specification of prior distributions; Section 0] provides a brief overview 
of an efficient computational algorithm for obtaining samples from the joint posterior, its 
implementation and methods for comparing goodness-of-fit across model specifications. 
Section [5] presents a comprehensive simulation study investigating the performa nce of the 


proposed framework, including a comparison with the methods of [Liquet et ah 


(2012). Sec¬ 


tion [6] reports on a detailed analysis of the motiving pancreatic cancer study; sensitivity 
analyses regarding the specification of certain model parameters are reported in Section 
0 Finally Section [S] concludes the paper with a discussion. Where appropriate, detailed 
derivations and additional results are provided in an online Supplementary Materials doc¬ 
ument. 


2 Risk of Readmission Among Patients Diagnosed with 
Pancreatic Cancer 


As mentioned in the 


2012 


CMS 


ntroduction, readmission is a key marker of quality-of-care flAsh et ah 


2013a.bl). To-date, however, studies of readmission have focused on health 


conditions that have relatively good prog nosis and/or low morta 


ure, m yocard ial infa rction and pneumonia ( Krumholz et al 


Epstein et ah 


, 2011 ). 


1997 


ity in cluding hear 


2011 


Jovnt et ah 


fail- 


2011 


Beyond these conditions, however, little is known about variation 
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in risk of readmission for patients diagnosed with terminal conditions such as pancreatic 
cancer. We are therefore currently engaged in a collaboration investigating readmission 
among Medicare enrollees diagnosed with pancreatic cancer. The overarching goals of the 
study are to improve end-of-life quality of care for these patients by first understanding 
patient-level risk factors associated with readmission and second understanding variation in 
risk at the level of the hospital (i.e. that not explained by differences in patient case-mix). 
Towards this we identified all n=5,298 Medicare enrollees who were diagnosed with pancre¬ 
atic cancer during a hospitalization at one of J=112 hospitals in the six New England states 
(Connecticut, Maine, Massachusetts, New Hampshire, Rhode Island and Vermont) between 
2000-2009. Information on the initial hospitalization and diagnosis, patient characteristics 
and co-morbid conditions, discharge destination and subsequent readmissions is obtained 
from the Medicare Fee-For-Service inpatient claims hie (Part A). Specific covariates of in¬ 
terest include sex (0/1 = male/fem ale), age, race ( 0/1 = white/non-white), the patients 
Charlson/Deyo comorbidity score (ISharab iani et a h. 2012), information on entry route for 
the initial admission (0/1 = from the ER/transfer from some other facility), whether or not 
the patient underwent a pancreatic cancer-specific procedure (resection, bypass, or stent), 
the length of hospitalization and the discharge destination. For the latter, patients could 
have been discharged to their home, their home with care, a hospice, an intermediate care or 
skilled nursing facility (ICF/SNF) or some other facility (e.g. a rehabilitation facility or to 
inpatient care). Table |T] provides a summary of observed distributions for these covariates. 

Also provided in Table [I] is a summary of the observed outcome information at 30 
and 90 days post-discharge. Specifically, each patient is classified into one of four groups: 
(1) they experienced a readmission event and were subsequently observed to die; (2) they 
experienced a readmission event but were censored prior to death; (3) they were observed 
to die without having experienced a readmission event; and, (4) they were censored prior to 
experiencing either a readmission or death event. The administrative censoring at 30 and 
90 days is driven by a several important factors. First, scientific and health policy interest 


regarding readmission has g enera. 


months following discharge flCMS 


ly focused on a patient’s experience in the immediate 


2013a]). The primary rationale for this is that post¬ 


discharge management for patients diagnosed with pancreatic cancer generally focuses on 
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Table 1: Covariate and outcome information for n=5,298 Medicare beneficiaries diagnosed with 
pancreatic cancer in the six New England states during a hospitalization between 2000-2009. 
Outcome information is considered with administrative censoring applied at 30 and 90 days post¬ 
discar ge. 




n 

Percent 

Covariate information 

Sex 

Female 

3,037 

57.3 


Male 

2,261 

42.7 

Age, years 

65-69 

727 

13.7 


70-74 

1,052 

19.9 


75-79 

1,226 

23.1 


80-84 

1,129 

21.3 


> 85 

1,164 

20.0 

Race 

White 

4,982 

94.0 


Non-white 

316 

6.0 

Charlson/Deyo comorbidity score 

< 1 

4,854 

91.6 


> 1 

444 

8.4 

Entry route 

Emergency room 

2,255 

42.6 


Transfer from another facility 

3,043 

57.4 

Procedure during hospitalization 

Yes 

1,291 

24.4 


No 

4,007 

75.6 

Length of hospitalization, days 

1-7 

3,170 

59.8 


8-14 

1,465 

27.7 


> 15 

663 

12.5 

Discharge destination 

Home 

1,823 

34.4 


Home with care 

1,571 

29.7 


Hospice 

419 

7.9 


SNF/ICF 

1,219 

23.0 


Other facility 

266 

5.0 

Outcome information with administrative censoring at 30 days 

Readmission and death 

205 

3.9 

Readmission and censored prior to death 


853 

16.1 

Death without readmission 


1,257 

23.7 

Censored prior to readmission or death 


2,983 

56.3 

Outcome information with administrative censoring at 90 days 

Readmission and death 

608 

11.5 

Readmission and censored prior to death 


930 

17.6 

Death without readmission 


1,912 

36.1 

Censored prior to readmission or death 


1,848 

34.9 


palliative care, with a specific emphasis on pain management. As patients and their health 
care providers coordinate this care, the early phases are particularly important for long- 





term success and are therefore of key interest. A second consideration is that readmission 
events that occur soon after a patient is discharged are more likely to be directly related 
to their diagnosis and subsequent care. Readmission events that occur a long time after 
diagnosis are less likely to be directly related to the quality of care they receive in the 
immediate aftermath of the diagnosis and, arguably, should not count against a hospitals 
performance. 


o 



Figure 1: Hospital-specific distributions of the outcome information for n=5,298 Medicare benefi¬ 
ciaries diagnosed with pancreatic cancer at one of J=112 hospitals in the six New England states 
between 2000-2009. Outcomes have been administratively censored at 90 days and distributions 
ordered according to the observed percentage of deaths (with and without readmission). 


A central feature of the Medicare data is that the n=5,298 patients are clustered within 
J=112 hospitals; cluster sizes vary from 10-420 with a median of 30 patients. The inherent 
clustering of patients within hospitals is important from both a statistical and a scientific 
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perspective: valid inference requires acknowledging potential correlation among patients 
and understanding between-hospital variation in readmission rates is a key scientific goal. 
Towards the latter, Figure Q] provides a barplot of the hospital-specific distributions of the 
four outcome groups based on censoring at 90 days. While there are many ways in which the 
<7=112 hospitals could be ordered, Figured] orders them according to the total percentage of 
patients readmitted within 90 days (i.e. with or without a subsequent death event). From 
the figure we see that there is substantial variation in observed readmission rates across 
hospitals, with the lowest being 5.6% and the highest being 64.3%. Moving beyond these 
raw adjusted rates would need to first account for case-mix differences across the hospitals, 
second account for death as a competing risk and third account for the cluster-correlation. 


3 A Bayesian Framework for Cluster-Correlated Semi- 
competing Risks Data 

3.1 Model specification 

Viewing ‘discharge’, ‘readmission’ and ‘death’ as three states, the underlying data gener¬ 


ating mechanism that gave rise to t 


rese data can b e represented by a multi-s 


speci fically an illness-death model (jAndersen et al 


1993 


Putter et al 


;ate model, 


2007 


Xu et. al, 


2010|). Letting Tj denote the time to non-terminal event and T 2 the time to the terminal 
event, the illness-death model is characterized by three hazard functions that govern the 
rates at which patients transitions between the states: 


hi (h) 

= lim-^Prpi e + A)| Ti 

A-U) A 

> t],T 2 > ti], 

O 

A 

-+o 

(1) 

h 2 (t 2 ) 

= bm — Pr[T 2 G [t 2 ,t 2 + A)\ Tj 
A-K) A 

> ^ 2 1 t 2 > t 2 ], 

t 2 > 0 

(2) 

^■3%2 \tl) 

= lim -J- Pr[T 2 e [t 2 , t 2 + A) | 7\ 
A-U) A 

— t\,T 2 > t 2 ], 

t 2 > t\ . 

(3) 


In practice, analyses based on the illness-death model characterized by (JTj)-(|3j) proceeds 
by placing structure on these functions, specifically as a function of covariates and frail¬ 
ties/random effects. Towards this, let Tjn and T ji2 denote the time to the non-terminal 
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event and time to the terminal event for the i th patient in the j th cluster, respectively, 
for i = 1 ,,rij and j = 1 ,,J. Furthermore, let X gig be a vector of time-invariant 
covariates for the i th patient in the j th cluster that will be considered in the model for the 
g th transition, g=l,2,3. Consider the following general modeling specification: 


r fjii X-jil') Vj 1 ) 

= 7 ji hmitjn) expjX^di + V)i}, 

tjn > 0 

(4) 

^ 2 (^ 2 ; Ijh Xji2, Vj2) 

= 7 ji h Q2 (t ji2 ) exp{Xj 2 /3 2 + V j2 }, 

tji 2 > 0 

(5) 

^ 3(^22 'Yjii jiS 1 ^j3) 

= la hwitji^tjn) exp {Xj a !3 3 + V j3 }, 

tji2 tjili 

( 6 ) 


where 7 is a shared patient-specific frailty, Vj = (Vji, Vj 2 , V) 3 ) is a vector of cluster- 
specific random effects, each specific to one of the three possible transitions, and (3 g is 
a transition -specific yector of fixed-effect log-hazard ratio regression parameters. As de¬ 


scribed by 


Xu ct ah 


(1201011 . model (Ed]) is often simplified in practice by either assuming 


that h 0 3 (tji 2 \tjn) = h 03 (tji 2 ) or that h Q3 (tj i2 \tjii) = h 03 (t ji2 -tjn). Given the former specifi¬ 
cation, the model is referred to as being Markov in the sense that the hazard for death given 
readmission does not depend on the actual time of readmission; under the latter specifica¬ 
tion, the model is referred to as semi-Markov. For simplicity we focus the exposition in this 
section on Markov models although note that the methods and computational algorithms 
have also been developed and implemented for the semi-Markov model; the analyses in 
Sections | 6 ] and [7] also consider both models. 


3.2 The observed data likelihood 


To complete the notation developed so far, let C gi denote the right censoring time for the 
2 th patient in the j th cluster. Furthermore, let Yjn = mm.(Tjn, T ]l2 , Cji ), A^i = 1 if Yjn 
= Tjn (i.e. a readmission event is observed) and 0 otherwise, Yj i2 = min (Tj i2 , Cji ) and, 
= 1 if Yji 2 = T ji2 (i.e. a death event is observed) and 0 otherwise. Finally, let Vji = 
{yjii,Sja,yji 2 ,5ji 2 } denote the observed outcome data for the 2 th patient in the j th cluster 
and H 0g (-) the cumulative baseline hazard function correspondin g to hn g (-). Let 7 and V 


denote the collections of the j gi and Vj, respectively. Following Putter et al. (120071 ). for a 
given specification of (|4])-((6]), the observed data likelihood as a function of the unknown 
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parameters $ = {/3i, (3 2 , fa, h 01 , h 02 , h 03 , 7 , V}, is given by: 


J n i 

apw = nn C(Vji\(3i, {3 2i hoi , ho2, ho3i 7^, Vj) 

j= 1 i=l 
J Tlj 

= n II { r rji h oi(yjii)vjn} Sjil(1 ~ 5ji2) {rfihoi(yjii)vjiih 0 3(yji2)vjiz} S]ll5jl2 
j= 1»=i 

X {7ji^02(2/7*2)^7722} (1 " <5jil)<5ji2 exp {-7 jiriyjn, y ji2 )} , ( 7 ) 

where rj jig = exp {Xj ig ( 3 g + V jg } and 2) = [i? 0 i(^1)77^1+ if 02 (ijii)^i2+ {-f 7 03 (^2) — 

77o3(tj ? :i)} 7 . 73 ] • 

In the remainder of this section, we complete the specification of the Bayesian model 
by providing detail on a range of possible choices for specification of the baseline hazard 
functions in (Rj)-(f6|), the population distribution for the hospital-specific random effects 
and, finally, prior distributions. To facilitate the exposition, Table [2] provides a summary 
of four possible specifications of the model along with the hyperparameters that require 
specification by the analyst. 

Table 2: Summary of four models based on parametric and non-parametric specifications of the 
baseline hazard functions and hospital-specific random effects distributions. Hyperparameters 
that require specification by the analyst are provided in parenthesis. Note, {ae,be), the hyper¬ 
parameters for the patient-specific frailty variance component, require specification in all of four 
models. 



Hospital-specific random effects, Vj 


MVN 

DPMt 

Baseline hazard functions, h 0g (-) 

{^v,Pv) 

(Vo,Po,t) 

Weibull 

Weibull-MVN 

Weibull-DPM 

PEM f (a K(i ,a^ g , b a , g ) 

PEM-MVN 

PEM-DPM 


t PEM: piecewise exponential model; DPM: Dirichlet process mixture 


3.3 Baseline hazard functions 

Within the frequentist paradigm estimation and inference for time-to-event models is often 
based on a partial likelihood which conditions on risk sets, removing the need for analysts 
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to specify baseline hazard functions. In the Bayesian paradigm, however, one is required to 
specify these functions. Here, we consider two strategies. The first assumes that the under¬ 
lying transition times follow Weibull(o! W)S , K w , g ) distributions, parameterized so that ho g (t) 
= a W) gK, W) gt aw ' 9 ~ l . While such a parametric specification is appealing due to its computa¬ 
tional simplicity, especially in small-sample settings, the Weibull is somewhat restrictive in 
that the corresponding hazard function is strictly monotone. As an alternative, we consider 
a non-parametric specification based on taking each of the log-baseline hazard function s 
to be a flexible mixture of piecewise constant functions flMcKeague and TighiouartL 2000). 
Briefly, let s 9imax denote the maximum observed time for transition g and partition (0, 
s 5 ,max] into K g + 1 intervals: s g = • • • , s g)Kg+ J, with s gfi = 0 and s g , Kg+1 = 


Jg, max* 


Given the partition (K g , s g ), we assume 


Ka +1 


log hog (t) = A g{t) = 1 b 9 , 




( 8 ) 


k=1 


where \ 9) k is the (constant) height of the log-baseline hazard function on the interval 
(s Si fc_i, s Si fc]. We refer to this specification as a piecewise exponential model (PEM) for 
the baseline hazard function. Note, w hile numerous options are available for specifying 


Ibrahim et al 


200l|), a key benefit of this structure is that it balances 


these functions (e.g. 

flexibility and computational convenience, since the integrals in the likelih ood ( spec ifically 
for the cumulative hazard functions) are replaced with summations flLee et all 120151 ). 


3.4 Hospital-specific random effects 

As with specification of the baseline hazard functions, we consider two options for the spec¬ 
ification of the population distribution for the J hospital-specific vectors of random effects. 
First, motivated by the current standard for analyses of readmission (i.e. a LN-GLMM), we 
consider a specification in which the Vj arise as i.i.d draws from a mean-zero multivariate 
Normal distribution with variance-covariance matrix Ey. The diagonal elements of the 
3x3 matrix Ey characterize variation across hospitals in risk for readmission, death and 
death following readmission, respectively, that is not explained by covariates included in 
the linear predictors. Crucially, that each random effect has its own variance component 
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allows the characterization of differential variation across hospitals for each of the three 
transitions. In addition, the off-diagonals of Ey permit covariation between the three ran¬ 
dom effects across the hospitals giving researchers the ability to characterize, for example, 
whether or not hospitals with high mortality rates tend to have low readmission rates. 

While conceptually simple and computationally convenient, the Normal distribution 
it is often criticized as being a strong assumption. As an alternative we consider the 
use of a so-called Bayesian nonparametric specification for the population distribution of 
Vj, spe cifically a Dirich i ef process mixture of multiva r iate Normal distributions (DPM- 
MVN) ( Ferguson , 


1973 


Bush an d MacEachernl, 


1996 


W alke r and Mallick , 


19971 ). One 


representation of this model is as follows: 


~ MVN 3 (/Xj, Ej), 

Hj , Ej | G ~ G, 

G ~ DP(G 0 ,r), (9) 


where Hj, E j are the cluster-specific latent mean and variance of Vj, which are taken to be 
draws from some (unknown) distribution G to which a Dirichlet process prior is assigned. 
Finally the Dirichlet process is indexed by Go, the so-called centering distribution, and r, 
the so-called precision parameter. 


3.5 Hyperparameters and prior distributions 

The proposed Bayesian framework is completed with the specification of prior distributions 
for unknown parameters introduced in Sections 13.1113.41 

3.5.1 Stage one parameters 

For each of the transition-specific regression parameters, /3 g for g= 1,2,3, a non-informative 
flat prior on the real line is adopted. For the shared patient-specific frailties, 7 ji, we assume 
that they arise from a Gamma(0" 1 , 6 1 " 1 ) distribution, parameterized so that if [ 7 ^] = 1 
and V[^fji] = 6. In the absence of prior knowledge on the frailty variance component, a 
Gamma(a 6 i, be) hyperprior for the precision 9~ l is adopted. 
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3.5.2 Baseline hazard functions 


For the parametric Weibull baseline hazard functions, since the hyperparameters have sup¬ 
port on (0, oo), we complete the specification of this model by adopting gamma prior distri¬ 
butions for both; that is, we take a w>g ~ Gamma(a a>g , b at9 ) and k W i9 ~ Gamma(a S)9 , b Kt9 ), 
g= 1,2,3. 

To complete the non-parametric PEM model specification, we specify that the K g + 1 
heights arise from a multivariate Normal distribution. Specifically, letting X g = .., 

A g: K g , ^ g ,K g +i) denote the transition-specific heights, we assume that \ g ~ MVN(/i Ag l, 
where /j,\ g is the overall mean, a\ g is a common variance component for the K g +1 
elements and S,\ 9 is a correlation matrix. To induce a priori smoothness in the baseline haz¬ 
ard functions we view the components of \ g in terms of a one-dimensional spatial problem, 
so that adjacent intervals can ‘borrow’ information from each other. To do t his we spec- 


ifv 5A via a Gaussian intrinsic conditional autoregression (ICAR) ([Besag and Koo perberg . 


19951 ). Additional technical details regarding the corresponding MVN-ICAR are provided 
in Supplementary Materials A. Finally, we specify a series of hyperpriors for the additional 
parameters introduced in the MVN-ICAR. In particular, we adopt a flat prior on the real 
line for p\ g and a conjugate Gamma(a CT)g , b a}9 ) distribution for the precision crjj 2 . For the 
ICAR specification, we avoid reliance on a fixed partition of the time scales by permitting 
the p artition (K g , s g ) to vary and be updated via a reversible jump MCMC scheme ( Green! . 
19951 ). Towards this we first adopt a Poisson( ck^) prior for the number of splits in the 
partition, K g . Conditional on the number of splits, we take the locations to be a priori 
distributed as the even-numbered order statistics: 


7T(s g \K g ) oc 


(2A-„+1)! nsr w - »m-i) 


( 10 ) 


Jointly, these choices form a time-homogeneous Poisson process prior for the partition (K g , 
s g ) so that a posteriori , after mixing over partitions as they arise in the MCMC scheme, the 
value of \ g (t) in any given small interval of time is characterized as a smooth exponentiated 


mixture of piecewise c onsta nt functions flArias and Gasbarral. 


2000 |; 


Haneuse et al 


20081) 


1994 


McKeague and Tighiouart 
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3.5.3 Hospital-specific random effects 

For the parametric specification of a single MVN 3 (0, 'Ey) distribution, we adopt a conju¬ 
gate inverse-Wishart (T,,, p v ) prior for the variance-covariance matrix 'Ey. Completion of 
the non-parametric DPM-MVN model requires specification of prior choices for the cen¬ 
tering distribution and the precision parameter. Here we take Go to be a multivariate 
Normal/inverse-Wishart (NIW) distribution for which the probability density function can 
be expressed as the product: 


/vxw(/ i ; Sl'ko, Po) — /mVAt(HO, £) X finv-Wishart^l^Q, P 0 ), 


where /dOI^d) is the density function for a distribution D indexed by G This cho ice is 
appe aling in that one can exploit prior-posterior conjugacy in the MCMC scheme flNeall . 


200 01). Finally, we treat hyperparameter the precision param eter in D PM-MVb 


tion, r, as unknown and assign a Gamma(a T , b T ) hyperprior flEscobar and West 


spec ifica- 


1995( 1 


4 Posterior Inference and Model Comparison 


4.1 Markov Chain Monte Carlo 

To perform estimation and inference for each of the models in Table [2] we use a random scan 
Gibbs sampling algorithm to generate samples from their joint posterior distributions. In 
the corresponding Markov chain Monte Carlo (MCMC) scheme, parameters are updated by 
either exploiting conjugacies inherent to the model structure or using a Metropolis-Hastings 
step. For models that adopt a PEM specification for the baseline hazard functions, updating 
the partition ( K g , s) requires a change in the dimension of the parameter space and a 


Metropolis-Hastings-Green step is used (IGree n, 


1995|). A detailed description of proposed 


computational scheme is given in Supplementary Materials B; as mentioned in Section IXT1 
the computation scheme has been developed for both the Markov and semi-Markov models 
for /z 03 (-). 

Finally , we note t hat the algorithms are implemented in the SemiCompRisks package for 


R (R Development Core Team! 20141. Given the complexity of the proposed models, and 
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the numerous updates in the MCMC scheme, C has been used as the primary computational 
engine to ensure that analyses can be conducted within a reasonable timeframe. 


4.2 Model comparison 

In practice, analysts have to balance model complexity with the realities of sample size 
and availability of information. While each of models in Table [2] has its own merit and 
utility, it may be of interest to directly compare their goodness of fit to the observed 
data. To this en d, we con si der two model assessment metrics: the deviance information 


criterion (DI C: ISpiegelhalter et al. 
tic (LPML; 


Geisser and E ddy. 


1979 


20021 ). and the log-pseudo marginal likelihood statis 


Gelfand and Mallickj, 


1995 ). Although DIC is often 


the default choice for model comparison in the Bayesian p aradi gm, its use in the context 
of complex hierarchical models requires care flCeleux et al.l . 120061 ). Specifically for models 
that condition on latent parameters, such as the patient-specific 7 ^ in models (J4J)-(J 6 J) , DIC 
computed on the basis of a likelihood that is marginalized with respect to these parameters 


performs more reliably as a metric 


likelihood that conditions on them (IM illar, 


or com pariso n than DIC computed on the basis of a 


2009]). For our purposes, since the Vy random 


effects are of intrinsic scientific interest, we propose to evaluate DIC and LPML on the 
basis of a partially marginalized likelihood, one that integrates solely over the distribution 
of the patient-specific frailties: 


J n i 


C*(V |<F*) = 


nn C(Vji |/3i, fa, An h 0 i, h 02 , h 0 3 , 7 , Vj)f e ( 7 ; e )dl (11) 

3 =1 i=l 


where <F* = {Ai, A2, A3, h 0 i, h 02l h 03 , 6,V}, C{'Dji\-) is given by expression (□) and /(•; 6) is 
the density of a Gamma(0 _1 , 0 _1 ) distribution (see Section [3.5.II) . 

Given expression (fill) , we therefore compute DIC as: 


DIC = D(<F*) + 2 p D , 


( 12 ) 


where /}($*) = (T) — 2 logZ/*('A,■^| < ^ , *) is the (marginal) deviance and <f>* is the posterior mean 




of <F*. The penalty term, p£>, is given by D($*) — £>($*), where D(<F*) is the posterior 
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mean of £)($*). Note, a model with smaller DIG indicates a better fit of the model for the 
data. 

The LPML statistic is computed as ^ log(CPO)jj, the sum of the logarithms of the 




patient-specific conditional predictive ordinate (CPO) flGeisserl . 119931 ). each defined as: 


CPOji = 


(13) 


where X)Gid denotes the data with the observation from the i th patient in the j th cluster 
removed. Intuitively, the CPO.,, is the posterior probability of the observed outcome for i th 
patient in the j th cluster, i.e. (yyi, 5yi, Uji 2 , dp/i), on the basis of a model fit to a dataset 
that excludes that particular patient. Thus, large values of CPO,,; attribute high posterior 
probability to the observed data and, therefore, indicate a better fit. Although a cl osed form 


expres sion for CPO.,, is not available for our proposed models, following 
(120001 ) we approximate CPOj, via a Monte Carlo estimator: 


Shao and Ibrahim 




(14) 


where {•h’d'h; q — 1,2,, Q} are MCMC samples drawn from the (marginal) joint posterior 
distribution of <f>*. 


5 Simulation Studies 


The performance of the proposed models is investigated through a series of simulation 
studies. The overarching goals of the simulation studies are to investigate the small sample 
operating characteristics of the models summarized in Table [2] un der a var i ety o f scenarios 


as well as to compare their performance with the methods of 


Liquet et al 


( 2012 ). 


5.1 Set-up and data generation 

Towards developing a comprehensive understanding of the performance of the proposed 
methods we consider six data scenarios that vary in terms of the true underlying baseline 
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hazard distributions, the true distribution of the cluster-specific random effects and the 
true extent of variation in the patient-specific frailties. Table E] provides a summary. In 
scenarios 1-5, the baseline hazard functions are set to correspond to the hazard of a Weibull 
distribution so that the event rates in the simulated data are similar to those in the ob¬ 
served Medicare data when the outcomes are administratively censored at t=90; specifically, 
we set (a Wt i, k 1U) i) = (0.8, 0.05), (a Wj2 , k w , 2 ) = (1.1, 0.01), and (a w , 3 , «„,,3)=(0.9, 0.01). To 
evaluate the performance of the model when the baseline hazard functions do not corre¬ 
spond to a Weibull, scenario 6 takes them to be piecewise linear functions: h 0g (t ) = {(k g - 
b g )t / 40+b g }I (£<40) + {(3kg-b g )/2-(kg-bg)t/80}I(t>40), with 6i=0.1, &2=0.05, 63 = 0 . 15 , and 
^1=^2=^3=0.0005 specified so that the true baseline hazard functions are not monotone 
increasing or decreasing functions like a Weibull. 


Table 3: Summary of six simulation scenarios explored in Section [3J 


Scenario 

Distribution of 
baseline hazard functions 

Distribution of 

cluster-specific random effects, V 3 

e 

1 

Weibull 

MVN(0, 0.25-7) 


0.50 

2 

Weibull 

/ 

4VN(0, 0.25-7) 
r 0.25 -0.10 -0.101 

\ 

1.00 

3 

Weibull 

MVN (0, 

-0.10 0.25 0.20 
-0.10 0.20 0.25 

) 

0.50 

4 

Weibull 

MVN(0, 0.25-7) 


0.00 

5 

Weibull 

0.5-MVN(0, 7)+0.5-MVN(0, 0.01-7) 

0.50 

6 

Piecewise linear 

MVN(0, 0.25-7) 


0.50 


With regard to the ‘true’ distribution of the cluster-specific random effects, scenarios 
1, 2, 4 and 6 consider a multivariate Normal distribution in which the components are 
independent. Scenario 3 expands on this by considering the impact of covariation across 
the Vj, while Scenario 5 examines the performance of the models when the true distribution 
is a mixture of two multivariate Normal distributions. 

Finally, with regard to the ‘true’ variance of the patient-specific frailties, scenarios 1, 
3, 5 and 6 consider a base value of 0=0.5. This value was chosen as a compromise across 
the posterior medians from the fits of the four models in Table [2] to the Medicare data (see 
Table |9] below). Scenario 2 considers the impact of greater variation in the patient-specific 
frailties, while Scenario 4 corresponds to a misspecification of the proposed model with the 
‘true’ 6 set to 0. 
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For each of the six scenarios we generated f?,=500 simulated datasets under the semi- 
Markov illness-death model described in Section 13.11 Across all simulated datasets, we 
set the number of clusters and cluster-specific sample sizes to be those observed in the 
Medicare data. Furthermore, we specified that each of the three transition-specific haz¬ 
ard functions depended on three covariates: Xj igl and X gig2 both Normal(0, 1) random 
variables and Xji 9} 3 a Bernoulli(0.5) random variable. The regression coefficients are set to 
/3 i=/32 =(0.5, 0.8, -0.5) and /33=(1.0, 1.0, -1.0), so that the covariate effects on the risk of 
the terminal event depend on whether or not the non-terminal event has occurred. Finally, 
we note that the function used to simulate the semi-competing risks data is available in 
the SemiCompRisks package. 

5.2 Analyses 

For each of the A=500 datasets under each of the six scenarios we fit each of the four models 
in Table [21 For the proposed models in which the baseline hazard function was specified via 
a Weibull distribution, we set (a Qj9 , b ag ) = (0.5, 0.01) and (a K , g , b K ^ g ) = (0.5, 0.05) for the 
transition-specific shape and rate parameters. For models in which a non-parametric PEM 
specification was adopted for the baseline hazard function, we set the prior Poisson rate 
on the number of intervals to be a g = 10. For the precision parameter in the MVN-ICAR 
specification, we set (a a , g , b aj9 ) = (0.7, 0.7) so that the induced prior for cr\ had a median 
of 1.72 and 95% central mass between 0.23 and 156. 

For the variance component associated with the patient-specific frailties, we set (a#, 
be ) = (0.7, 0.7); that is the same prior was used for the precision O^ 1 for the 7 ^ frailties 
as the precision component in the MVN-ICAR specification for the PEM model. For the 
hospital-specific random effects variance components, given a MVN specification, we set 
(\&„, p v ) = (/ 3 , 5) so that the induced prior on Ey has a prior mean given by the 3x3 
identity matrix. The same prior was adopted for the variance-covariance matrix of the 
centering distribution of the DPM-MVN specification, G 0 ; that is, we set (Tq, p 0 ) = (-G, 
5). Finally, for the precision parameter in the DPM-MVN specification we set (a T , b T ) = 
(1.5, 0.0125) so that a priori r had a mode of 40 and standard deviation of 98. Given 
the prior specifications, two independent chains were run for a total of six million scans 
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each; the Gelman-Rubin potential scale reduction (PSR) statistic ( Gelman et all 2013) 


was used to assess convergence, specifically requiring the PSR to be less than 1.05 for all 
model parameters. 

In additio n to t he m odels in Table El we analyzed each simulated dataset using the 


methods of 


Liquet et al 


(2012). Specifically, we con sid ered the ‘shared frailty’ (SF) model 


implemented in the f railtypack package for R ({Ro ndeau et ah, 2012) and summarized us¬ 
ing notation consistent with that adopted in this manuscript in Supplementary Materials 
Section C. Briefly, this model adopts a Cox-type regression structure for each transition- 
specific hazard, as we do in expressions (141) - (l6l) . For the baseline hazard functions, two 
options are available: one that corresponds to a Weibull distribution and another where 
each hg(-) is specified via a flexible penalized smoothing spline. To distinguish these mod¬ 
els, we refer to them as the Weibull-SF and Spline-SF models, respectively. In contrast 
to the specification in (ITl)-(l6j). the SF model introduces a cluster-specific frailties as a mul¬ 
tiplicative factors for each transition-specific hazard. Two options are available for the 
distribution of these factors across the clusters; either they arise from three independent 
gamma distributions or they arise from three independent log-Normal distributions. For 
either option, estimation and inference is performed within the frequentist paradigm specif¬ 
ically based on an integrated likelihood that marginalizes out the cluster-specific frailties; 
estimation of the latter is performed via empirical Bayes. In this paper, we present the 
results from the SF models that adopt independent gamma distributions for cluster-specific 
frailties while we provide those from the SF models with independent log-Normal frailties 
in Supplementary Materials D. Finally, we note that in contrast to the specification in 
expressions (jl])-©, the SF model does not account for within-patient correlation. That 
is there is no quantity that corresponds to the patient-specific jji terms in the proposed 
models. 


5.3 Results 

5.3.1 Baseline survivor functions 

Figure El presents the mean estimated transition-specific baseline survival functions under 
scenarios 1, 4 and 6 across the six models. Under scenarios 1 and 4, for which the baseline 
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hazard functions are Weibnll, all four of the proposed models estimate the three baseline 
survivor functions very well. In contrast the two SF models only perform well in scenario 
4 for which 6=0. This is to be expected since, as described in detail in Supplementary 
Materials Section C, the SF model does not include patient-specific frailties; effectively, 
it assumes that 6=0 even when it is not. In scenario 6, for which the baseline hazard 
functions are not Weibnll, the proposed PEM-MVN and PEM-DPM specifications capture 
the true shape of the baseline survivor functions well; all four of the models that assume 
the baseline hazard function to be a Weibull, however, are unable to capture the shape. 

5.3.2 Regression parameters and 6 

Focusing on scenarios 1-3, each corresponding to a ‘true’ Weibull-MVN model, Table [4] 
indicates that all four of the proposed models in Table [2] perform very well in terms of 
estimation and inference for / 3 \ and 6. Across the board, we find that percent bias is 
no larger than 3.2% and the estimated coverage probabilities are all close to the nominal 
0.95. In contrast, both the Weibull-SF and Spline-SF models yield point estimates of (3\ 
that are significantly biased and, as such, have poor coverage probabilities. The poor 
performance of the SF models is likely tied to the fact that they do not account for within- 
patient correlation; hence 6 is not estimated by these models. The results for these models, 
however, is dramatically improved under scenario 4 for which the true value of 6 is zero (i.e. 
the scenario they explicitly accommodate). Interestingly, the four proposed models each 
exhibit a small amount of bias under this scenario (up to approximately 5%). In addition 
the coverage probabilities for /3 n and f3\2 are poor, particularly for the two models that 
adopt a PEM specification for the baseline hazard function. In scenario 5, we again see 
that all four of the proposed models perform well. Finally, under scenario 6 we see that 
the PEM-MVN and PEM-DPM models perform very well in terms of bias and coverage. 
In contrast, the Weibull-MVN and Weibull-DPM models perform poorly, particularly with 
respect to estimation of 9, illustrating the potential danger of adopting a parametric Weibull 
baseline hazard function when the truth is not a Weibull. 

While Table [4] explores estimation and the (valid) quantification of uncertainty, Table 
[5] examines the relative merits of the various analysis approaches in terms of efficiency. 
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Table 4: Estimated percent bias and coverage probability for j3\ and 0 for six analyses described in Section 15.21 across six simulation 
scenarios given in Table [3l Throughout values are based on results from i?=500 simulated datasets. 


Percent Bias Coverage Probability 


Scenario 


True 

Weibull 

Weibull 

Weibull 

PEM 

PEM 

Spline 

Weibull 

Weibull 

Weibull 

PEM 

PEM 

Spline 



value 

-MVN 

-DPM 

-SF 

-MVN 

-DPM 

-SF 

-MVN 

-DPM 

-SF 

-MVN 

-DPM 

-SF 


Pu 

0.50 

0.1 

0.2 

-19.8 

0.4 

0.4 

-21.0 

0.96 

0.96 

0.01 

0.95 

0.96 

0.00 

1 

Hl2 

0.80 

0.2 

0.3 

-19.7 

0.5 

0.4 

-21.0 

0.95 

0.95 

0.00 

0.96 

0.97 

0.00 


fill 

-0.50 

0.3 

0.3 

-19.8 

0.3 

0.3 

-21.2 

0.97 

0.96 

0.31 

0.96 

0.96 

0.25 


0 

0.50 

1.0 

1.3 


1.4 

1.2 


0.95 

0.95 


0.93 

0.94 



fin 

0.50 

-0.1 

-0.0 

-31.8 

0.1 

0.1 

-32.8 

0.94 

0.94 

0.00 

0.94 

0.93 

0.00 

2 

Hl2 

0.80 

0.1 

0.2 

-31.7 

0.4 

0.3 

-32.7 

0.97 

0.97 

0.00 

0.94 

0.95 

0.00 


fill 

-0.50 

1.2 

1.3 

-31.1 

1.1 

1.1 

-32.2 

0.94 

0.95 

0.05 

0.94 

0.94 

0.04 


0 

1.00 

0.4 

0.7 


0.7 

0.6 


0.94 

0.95 


0.94 

0.95 



fin 

0.50 

0.3 

0.3 

-19.9 

0.7 

0.7 

-21.0 

0.94 

0.94 

0.00 

0.93 

0.94 

0.00 

3 

Hl2 

0.80 

0.4 

0.4 

-19.8 

0.8 

0.8 

-20.9 

0.94 

0.94 

0.00 

0.94 

0.94 

0.00 


Pii 

-0.50 

0.4 

0.3 

-20.1 

0.5 

0.6 

-21.2 

0.96 

0.96 

0.31 

0.95 

0.96 

0.27 


0 

0.50 

2.0 

2.1 


3.2 

3.2 


0.96 

0.96 


0.93 

0.95 



fin 

0.50 

3.7 

3.7 

0.2 

4.7 

4.6 

0.3 

0.87 

0.86 

0.96 

0.81 

0.83 

0.96 

4 

P 12 

0.80 

3.6 

3.6 

-0.0 

4.5 

4.5 

0.1 

0.80 

0.79 

0.95 

0.69 

0.70 

0.95 


Hl3 

-0.50 

4.0 

4.0 

0.2 

4.8 

4.7 

0.2 

0.93 

0.94 

0.94 

0.93 

0.93 

0.93 


0 

0.00 














Hi 1 

0.50 

-0.3 

0.1 

-20.3 

0.0 

0.3 

-21.1 

0.94 

0.95 

0.00 

0.96 

0.96 

0.00 

5 

Hl2 

0.80 

0.0 

0.3 

-20.0 

0.3 

0.6 

-20.9 

0.95 

0.95 

0.00 

0.96 

0.96 

0.00 


Hll 

-0.50 

-0.2 

0.2 

-20.4 

-0.2 

0.2 

-21.3 

0.94 

0.94 

0.29 

0.94 

0.94 

0.25 


0 

0.50 

-0.2 

1.0 


0.4 

1.3 


0.95 

0.95 


0.95 

0.96 



fill 

0.50 

9.3 

9.4 

-22.1 

0.4 

0.3 

-25.9 

0.58 

0.57 

0.00 

0.94 

0.94 

0.00 

6 

Hl2 

0.80 

9.7 

9.8 

-22.0 

0.5 

0.5 

-25.8 

0.20 

0.20 

0.00 

0.94 

0.95 

0.00 


Hii 

-0.50 

10.2 

10.2 

-21.6 

0.8 

0.7 

-26.1 

0.81 

0.80 

0.21 

0.93 

0.94 

0.10 


0 

0.50 

52.8 

53.0 


1.8 

1.7 


0.00 

0.00 


0.95 

0.96 













Baseline survival function Baseline survival function Baseline survival function 
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Figure 2: Estimated transition-specific baseline survival functions, So g (-) = exp(-Ho g (-)), for 
each six analyses described in Section 15.21 under simulation scenarios 1, 4 and 6. 
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Table 5: Average relative width of 95% credible/confidence intervals for /3\ and 9, with the 
Weibull-MVN model taken as the referent, across six simulation scenarios given in Table [3l 
Throughout values are based on results from i?=500 simulated datasets. 


Scenario 


Weibull 

Weibull 

Weibull 

PEM 

PEM 

Spline 



-MVN 

-DPM 

-SF 

-MVN 

-DPM 

-SF 


Pn 

1.00 

1.00 

0.81 

1.02 

1.02 

0.81 

1 

@12 

1.00 

1.00 

0.77 

1.04 

1.04 

0.77 


@13 

1.00 

1.00 

0.84 

1.00 

1.01 

0.83 


9 

1.00 

1.00 


1.10 

1.12 



@n 

1.00 

1.00 

0.73 

1.02 

1.02 

0.73 

2 

@12 

1.00 

1.00 

0.69 

1.03 

1.04 

0.69 


@13 

1.00 

1.00 

0.76 

1.00 

1.00 

0.76 


9 

1.00 

1.00 


1.12 

1.14 



@11 

1.00 

1.00 

0.81 

1.02 

1.02 

0.81 

3 

P 12 

1.00 

1.00 

0.76 

1.04 

1.04 

0.77 


@13 

1.00 

1.00 

0.83 

1.00 

1.01 

0.83 


9 

1.00 

1.00 


1.10 

1.13 



P U 

1.00 

1.00 

0.95 

1.02 

1.01 

0.96 

4 

@12 

1.00 

1.00 

0.94 

1.03 

1.03 

0.95 


@13 

1.00 

1.00 

0.96 

1.01 

1.01 

0.96 


9 

1.00 

1.00 


1.09 

1.09 



@11 

1.00 

1.00 

0.81 

1.02 

1.02 

0.81 

5 

@12 

1.00 

1.00 

0.77 

1.03 

1.03 

0.77 


@13 

1.00 

1.00 

0.83 

1.00 

1.00 

0.83 


9 

1.00 

1.00 


1.09 

1.09 



@11 

1.00 

1.00 

0.74 

0.94 

0.95 

0.73 

6 

@12 

1.00 

1.00 

0.72 

0.96 

0.97 

0.71 


@13 

1.00 

1.00 

0.76 

0.93 

0.93 

0.75 


9 

1.00 

1.00 


0.89 

0.90 



Specifically, we computed the average relative width of 95% credible/conhdence intervals 
for (3\ and 6 under each analysis with the Weibull-MVN model taken as a common ref¬ 
erent. Comparing the Weibull-DPM to the Weibull-MVN as well as the results between 
the PEM-MVN and PEM-DPM we see that there is no loss of efficiency for any of the 
regression parameters, and minimal loss for 6 , if one adopts the flexible DPM specification 
for the cluster-specific random effects, even if the true distribution is a MVN. Under all 
five scenarios for which the true baseline hazard functions were Weibull hazard functions, 
the two models that adopt a PEM specification have somewhat wider credible intervals 
particularly for 9. However, as expected, the 95% credible intervals for the two PEM mod¬ 
els under scenario 6 are somewhat tighter indicating improved efficiency when the true 
baseline hazard functions are not Weibull hazard functions. Finally, across all scenarios, 
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the estimated 95% confidence intervals for the two SF models are substantially tighter than 
those for any of the proposed analyses, although this must be balanced with the high bias 
shown in Table |4j 


5.3.3 Cluster-specific random effects 

Finally, Table [H] investigates the relative performance of the various analyses with respect 
to estimation of the cluster-specific random effects. Specifically, we calculated the mean 
squared error of prediction (MSEP) given by: 


1 

RJ 


R J 


r =1 j =1 



(15) 


where V r j g is the cluster-specific random effect for the j th cluster in the transition g for the 
r th simulated data set, r= 1,... ,R. For each of the four proposed models, V r j g was taken as 
the corresponding posterior median. For the two SF models, V r j g was taken as a the log of 
the empirical Bayes estimates of the transition/cluster-specific frailties (see Supplementary 
Materials C for details). We note, however, that for some of the simulated datasets, the 
empirical Bayes estimates returned by the current implemented in the f railtypack package 
were zero. Since taking the log of these estimates would yield V r j g = —oo, we calculated 
MSEP over the random effects for which the empirical Bayes estimate was non-zero; to 
place these values in context, Table |H] also reports the percentage of instances where a 
frailty was estimated to be zero. 

From Table [6] we see that under scenarios 1-4, for which the true model is a Weibull-MVN 
model, the Weibull-MVN analysis generally performs the best. Comparing the Weibull- 
MVN and PEM-MVN results across these scenarios, we see that over-specification of the 
baseline hazard functions (i.e. adoption of the more flexible PEM specification) does not 
meaningfully impact MSEP. In addition, comparing the Weibull-MVN and Weibull-DPM 
results we see that over-specification of the random effects structure (i.e. adoption of the 
more flexible DPM specification) does not adversely affect MSEP either. When the true 
distribution of the random effects is not a multivariate Normal distribution, however, as in 
scenario 5, both the Weibull-DPM and PEM-DPM models outperform their MVN coun- 
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Table 6: Mean squared error of prediction (xl0~ 2 ) for cluster-specific random effects based on six 
analyses described in Section 15.21 across six data scenarios given in Table [3j Throughout values 
are based on results from R =500 simulated datasets. 


Scenario 


Weibull 

-MVN 

Weibull 

-DPM 

Weibull 

-SF 

%Ft 

PEM 

-MVN 

PEM 

-DPM 

Spline 

-SF 

%F 


Pi 

5.25 

5.27 

6.40 


5.27 

5.27 

6.39 


1 

e 2 

7.66 

7.70 

8.70 

17.8 

7.67 

7.72 

8.68 

0.2 


e 3 

9.91 

9.95 

12.13 


9.91 

9.96 

12.11 



Ei 

6.36 

6.41 

8.10 


6.37 

6.41 

8.09 


2 

e 2 

8.76 

8.85 

10.23 

10.4 

8.77 

8.86 

10.20 

0.0 


e 3 

11.13 

11.19 

13.85 


11.13 

11.19 

13.91 



E x 

5.03 

5.04 

6.27 


5.04 

5.04 

6.22 


3 

e 2 

6.34 

6.34 

8.28 

15.8 

6.36 

6.36 

8.24 

0.0 


e 3 

7.55 

7.49 

11.66 


7.57 

7.55 

11.69 



Ei 

3.84 

3.85 

4.99 


3.87 

3.87 

5.01 


4 

e 2 

6.25 

6.27 

7.19 

12.8 

6.25 

6.27 

7.12 

0.4 


e 3 

7.89 

7.90 

9.57 


7.90 

7.91 

9.52 



Ei 

6.95 

6.26 

10.87 


6.96 

6.27 

10.86 


5 

e 2 

11.52 

10.50 

14.95 

12.8 

11.50 

10.52 

14.92 

0.2 


e 3 

15.46 

14.66 

25.04 


15.46 

14.72 

24.94 



E x 

5.05 

5.01 

6.34 


4.89 

4.85 

6.26 


6 

e 2 

7.58 

7.55 

8.60 

5.4 

7.41 

7.39 

8.49 

1.4 


e 3 

6.72 

6.65 

13.42 


6.44 

6.40 

13.70 



t % of times SF models yield at least one of Vj being — oo, resulting in MSEP being oo 


terparts, illustrating the potential benefit of the more flexible DPM specification. Further¬ 
more, when the true baseline hazard functions do not correspond to a Weibull distribution, 
the MSEP for the two PEM models are, as expected, smaller than the corresponding val¬ 
ues for the two Weibull models, illustrating the potential benefit of the more flexible PEM 
specification. Finally, we find that the empirical Bayes estimates of the cluster-specific 
random effects from the SF models perform relatively poorly when compared to the corre¬ 
sponding estimates from the proposed methods. For example, the Splinc-SF model yields 
approximately 14% to 55% higher MSEP than our proposed PEM-MVN model across the 
six scenarios. 
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6 Analysis of Medicare Data 


6.1 Analysis details and prior specifications 

Returning to the motivating application of readmissions following a diagnosis of pancre¬ 
atic cancer, we fit each of the four models summarized in Table [2] to the Medicare data 
under both the Markov and semi-Markov assumption for hs(-) (see Section I3T| . Based on 
the rationale provided in Section [21 we administratively censored observation time at 90 
day. Given the results from the simulation studies, specifically with respect to estimation 
of the cluster -speci fic random effects, we decided not to fit the shared frailty models of 


Liq uet et ah 


( 2012 ). We did, however, perform an analysis based on a LN-GLMM model 


since this model is the current standard for analyzing variation in the risk of readmission 
and we believed it would be instructive to examine the potential impact of ignoring death 
as a competing force. Towards this, let Y* = 0/1 be a binary indicator of whether or not 
the i th patient in the j th hospital readmitted within 90 days of discharge. Note, if a patient 
died prior to readmission within 90 days their outcome was set to Y* = 0. The LN-GLMM 
is then given by: 


logitPr(K,- = l|V-„U) = Xgp + V, 


(16) 


where V* is a hospital-specific random effect for readmission taken to be Normally dis¬ 
tributed with mean zero and a constant variance, a 2 . To complete the Bayesian specifica¬ 
tion of this model, we adopted a Gamma(0.7, 0.7) prior for the precision <j~ 2 . For the four 
proposed models, the hyperparameters outlined in Table [2] are specified as in Section 15.21 
Throughout the analyses, to ensure the baseline hazard functions in the proposed mod¬ 
els and the (overall) intercept in the LN-GLMM retained reasonable interpretations, age 
was standardized so that ‘zero’ corresponded to age 77 years and a one-unit increment 
corresponded to a 10-year contrast. Furthermore, length of stay during the initial hos¬ 
pitalization was also standardized so that ‘zero’ corresponded to 10 days and a one-unit 
increment corresponded to a 7-day contrast. 
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6.2 MCMC 


Towards obtaining summaries of the joint posterior distributions we ran 3 independent 
chains of the proposed MCMC scheme, each for a total of 6 million scans. Convergence 
was evaluated by inspection of trace plots as well as calculation of the PSR statistic; an 
MCMC scheme was determined to have converged if the PSR statistic was less than 1.05 for 
all parameters in the model (see Supplementary Materials E). Although the hierarchical 
models are complex and include a large number of parameters, the proposed algorithm 
achieved an overall acceptance rate of 35% across the various Metropolis-Hastings and 
Metropolis-Hastings-Green steps. To provide a sense of computational time, the most 
complex of our proposed models (i.e. the PEM-DPM model), the implementation in our 
E package is able to generate 1 million scans in 30 minutes on a 2.5 GHz Intel Core i7 
MacBook Pro; for the least complex of the models (i.e. the Weibull-MVN model), the 
implementation is able to generate 1 million scans in 10 minutes on the same machine. 


6.3 Results 


6.3.1 Overall model fit 

Table [7] provides DIC and LPML for the eight model fits. For the DIC measure, a general 
rule of thumb for model comparison is to consider differences of less than 2 to be negligible, 
differences between 2 and 6 to indicative of positive support for the model with the lower 
value and differences greater than 6 to be stron g support in value of the model with the lower 


value flSpiegelhalter et all 


2002 


Mi llarl . 20Q9J). For LPML, one can compute the so-called 


pseudo Bayes 


actor (PBF) for two models by exponentia ting difference in their LPML 


values (j Hanson . 200.6). While the conventional Bayes factor I Kass and Raftem 119951 1 tends 


to find which model explains the observed data best, predictive methods such as PBF 
attempt to find which model gives the best predictions for future ob servati ons when the 
same process as the original data is used to generate the observations flKadane and Lazar . 
20041 1 

Based on pairwise comparisons of the values in Table |7] we draw a number of conclu¬ 
sions. First, each of the models in which a semi-Markov specification is made for hos(-) 
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Table 7: DIC and LPML for eight models fit to the New England Medicare data. 



DIC 

LPML 

Weibull-MVN 

46184.3 

-23101.6 

Markov Weibull-DPM 

46174.1 

-23101.2 

PEM-MVN 

45609.2 

-22812.6 

PEM-DPM 

45606.8 

-22810.7 

Weibull-MVN 

46163.7 

-23088.8 

semi-Markov Weibull-DPM 

46153.0 

-23086.7 

PEM-MVN 

45574.1 

-22790.9 

PEM-DPM 

45569.0 

-22789.3 


has a substantially better fit to the data than the corresponding model in which a Markov 
assumption is made for h 03 (-); differences in DIC and the PBF range between 20.6-37.8 
and the order of 10 5 -10 9 , respectively. Second, both DIC and LPML indicate that models 
for which a PEM specification was adopted for the baseline hazard functions have substan¬ 
tially better fit to the data than models for which a Weibull hazard function was adopted; 
differences in DIC and the PBF range between 567.3-589.6 and the order of 10 125 -10 129 , 
respectively. Finally, although DIC indicates a somewhat better fit for models that adopt 
a DPM for the random effects distribution compared to a MVN specification, the LPML 
values are less convincing in this regard; differences in DIC and the PBF range between 
2.4-10.7 and 1.5-8.2, respectively. 

6.3.2 Baseline survival functions 

Since hazard functions are notoriously difficult to interpret, Figure [3] provides estimates 
of the corresponding baseline survival functions. Specifically, they provide pointwise time- 
specific posterior medians for Sq 9 (■) for a 77-year old white female patient who had a 
Charlson/Deyo comorbidity score of 0 or 1, whose initial hospitalization lasted 10 days 
and during which they had no pancreatic cancer-related procedures, and were eventually 
discharged to their own home. In panels (a)-(c) results are presented for models for which 
a Markov assumption was adopted for hos(-); panels (d)-(f) present results for models for 
which a semi-Markov assumption was adopted. 

From panels (a) and (d) we see that all eight models indicate similar risk of readmis¬ 
sion within the first 30 days. After 30 days, however, models with a parametric Weibull 
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Readmission 


Death without readmission 


Death after readmission 


(a) 



(b) 



(c) 



0 30 60 90 


30 60 90 


30 60 90 


Time since discharge, days 


Time since discharge, days 


Time since discharge, days 



Time since discharge, days 


Time since discharge, days 


Time since readmission, days 


Figure 3: Estimated transition-specific baseline survival functions, So g (-) = exp(— Ho g (-)), for the 
eight models fit to the New England Medicare data (see Table CD- Panels (a)-(c) are estimates 
based on a Markov specification for /io3(■); panels (d)-(f) are estimates based on a semi-Markov 
specification. Note the difference in time scale between panels (c) and (d) due the differing 
specifications for hos(-). 


specification for the baseline hazard functions indicate substantially higher overall risk for 
readmission. Note, most of the observed readmission events occur relatively soon after 
discharge with a median of 18 days and 75% of observed events occurring within 40 days. 
As such, the posterior mass is being assigned to values of the two Weibull hyperparameters, 
(a g , K g ), that fit the early time periods well to the detriment of fitting late periods relatively 
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poorly. From panels (b) and (e) a similar phenomenon is observed for the baseline survival 
function for death without readmission for which the median event time is 20 days and, 
again, approximately 75% of observed events occurring within 40 days. In contrast, since 
the distribution of time to death following readmission is more spread out (median=43 
days, IQR=40 days) the estimated baseline survival functions under the Weibull and PEM 
are more similar (see panels (c) and (f)). 

6.3.3 Regression parameters 

Posterior summaries for the vector of hazard ratio (HR) parameters for readmission, exp(/3i), 
are presented in Table [HI For brevity, based in part on the conclusions drawn from Table 

[7] results are only presented for models for which a semi-Markov specification was adopted 
for h 03 (-); additional results, particularly for exp(/3 2 ) and exp(/3 3 ) are provided in Supple¬ 
mentary Materials E. In addition, posterior summaries for the vector of odds ratio (OR) 
parameters, exp(/3*) from model flTTjh are also presented. 

Recognizing that the interpretation of the HR and OR parameters differ (due to the 
different set of frailties/random effects that are conditioned upon), the results in Table 

[8] indicate that the LN-GLMM qualitatively identifies a different set of risk factors for 
readmission than the results based on the proposed framework. For instance, while there 
is evidence of lower risk for readmission among females diagnosed with pancreatic cancer 
under the semi-competing risks approach (e.g. HR 0.80; 95% Cl 0.70, 0.90 in Weibull- 
MVN), one cannot draw the same conclusion based on the LN-GLMM (OR 0.91; 95% Cl 
0.80, 1.03). In addition, under the LN-GLMM model there is no evidence of a relationship 
between source of entry to the initial hospitalization (OR 0.99; 95% Cl 0.86, 1.14) while 
under each of the semi-competing risks analysis models there is evidence that patients 
who enter the initial hospitalization via some route other than the emergency room are at 
higher risk of readmission (e.g. HR 1.12; 95% Cl 1.00, 1.28). Conflicting results are also 
found with respect to discharge destination. In particular, under the LN-GLMM model 
patients who are discharged to home with care, a hospice, a ICF/SNF or some ‘other’ 
facility (e.g. a rehabilitation center) have statistically significant lower estimates odds of 
readmission than patients discharged to their home without care. In contrast, results from 
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Table 8: Posterior summaries (medians and 95% credible intervals (Cl)) for hazard ratio (HR) parameters for readmission, exp(/3i), 
from semi-competing risks data analyses and odds ratios (OR) based on a LN-GLMM. For the latter, results are reported for models 
for which a semi-Markov specification was adopted for /io3 (') 



LN-GLMM 

OR (95% Cl) 

Weibull-MVN 
HR (95% Cl) 

Weibull-DPM 
HR (95% Cl) 

PEM-MVN 

HR (95% Cl) 

PEM-DPM 

HR (95% Cl) 

Sex 






Male 

1.00 

1.00 

1.00 

1.00 

1.00 

Female 

0.91 (0.80, 1.03) 

0.80 (0.70, 0.90) 

0.79 (0.70, 0.90) 

0.85 (0.76, 0.96) 

0.85 (0.76, 0.95) 

Age! 

0.87 (0.83, 0.92) 

0.90 (0.86, 0.95) 

0.90 (0.86, 0.95) 

0.91 (0.87, 0.94) 

0.91 (0.87, 0.95) 

Race 






White 

1.00 

1.00 

1.00 

1.00 

1.00 

Non-white 

1.17 (0.90, 1.51) 

1.12 (0.86, 1.45) 

1.12 (0.86, 1.46) 

1.11 (0.89, 1.40) 

1.12 (0.89, 1.38) 

Source of entry to initial hospitalization 






Emergency room 

1.00 

1.00 

1.00 

1.00 

1.00 

Other facility 

0.99 (0.86, 1.14) 

1.18 (1.03, 1.36) 

1.19 (1.03, 1.35) 

1.12 (1.00, 1.27) 

1.13 (1.00, 1.28) 

Charlson/Deyo score 






< 1 

1.00 

1.00 

1.00 

1.00 

1.00 

> 1 

1.31 (1.05, 1.64) 

1.50 (1.19, 1.85) 

1.49 (1.20, 1.85) 

1.41 (1.15, 1.68) 

1.40 (1.15, 1.70) 

Procedure during hospitalization 






No 

1.00 

1.00 

1.00 

1.00 

1.00 

Yes 

0.73 (0.61, 0.86) 

0.45 (0.38, 0.53) 

0.44 (0.37, 0.54) 

0.56 (0.48, 0.66) 

0.57 (0.48, 0.66) 

Length of stay* 

1.14 (1.06, 1.22) 

1.15 (1.07, 1.24) 

1.15 (1.07, 1.24) 

1.12 (1.05, 1.18) 

1.12 (1.05, 1.19) 

Discharge location 






Home without care 

1.00 

1.00 

1.00 

1.00 

1.00 

Home with care 

0.68 (0.58, 0.79) 

0.95 (0.82, 1.11) 

0.96 (0.82, 1.11) 

0.90 (0.79, 1.02) 

0.89 (0.78, 1.02) 

Hospice 

0.06 (0.04, 0.10) 

0.39 (0.22, 0.62) 

0.39 (0.22, 0.62) 

0.27 (0.15, 0.45) 

0.27 (0.16, 0.43) 

ICF/SNF 

0.44 (0.36, 0.53) 

0.88 (0.72, 1.07) 

0.88 (0.73, 1.08) 

0.76 (0.64, 0.90) 

0.76 (0.65, 0.90) 

Other 

0.56 (0.41, 0.76) 

1.04 (0.75, 1.42) 

1.06 (0.75, 1.44) 

0.91 (0.70, 1.18) 

0.90 (0.68, 1.18) 


t Standardized so that 0 corresponds to an age of 77 years and so that a one unit increment corresponds to 10 years 
* Standardized so that 0 corresponds to 10 days and so that a one unit increment corresponds to 7 day 





the semi-competing risks analyses fail to indicate differences between patients discharged 
to home without care and those discharged to home with care (e.g. ffR 0.90; 95% Cl 
0.79, 1.02) or to some other facility. Furthermore, while patients discharged to either a 
hospice or ICF/SNF have significantly lower odds of readmission, the estimated effects 
are substantially attenuated (e.g. compare OR=0.06 under the LN-GLMM to ffR=0.27 
under the PEM-MVN model). Finally, consistent with the assessment of model fit in Table 
[7l Table [ 8 ] indicates that for estimation and inference for regression parameters differs 
somewhat between models based on a Weibull baseline hazard specification and models 
based on a PEM specification. Comparing the Weibull-MVN model to the PEM-MVN 
model, for example, estimates for gender, Charlson/Deyo score and whether or not the 
patient underwent a procedure during the hospitalization are all attenuated; in contrast 
estimates for discharge location are generally strengthened under the PEM-MVN model, 
in some cases achieving statistical significance. 


6.3.4 Variance components 


Table [9] provides posterior summaries for the standard deviation of the patient-specific 
frailty distribution, yf§, as well as components of the variance-covariance matrix for the 
hospital-specific V = {V\ , V 2 , V 3 ) from models in which a semi-Markov specification was 
adopted for h 03 (-). For the latter, the summaries are directly with respect to the components 
of £y under a MVN specification; under the two DPM specifications, posterior summaries 
are reported for the marginal total variance-covariance matrix obtained by applying the 
law of total cu mulan ce: )T)Li{(/ to, - p) (/to,, - p) T + £where p = Y^j=i PmJJ 


(lOhlssen et al. 


20071 ). From the Table we see that the components of variation (particu¬ 


larly the standard deviation components) are generally smaller in magnitude for models in 
which a PEM specification for the baseline hazard functions was adopted. For example, un¬ 
der the Weibull-MVN model the posterior median of \/6 is 1.03, whereas the corresponding 
posterior median under the PEM-MVN model is 0.61. This is likely due to the 7 ^ patient- 
specific frailties not only representing patient-level heterogeneity but also accounting, in 
part, for misspecihcation of the Weibull model when the underlying baseline hazard func¬ 
tions are not Weibull. Qualitatively, across all four model specifications, we find that there 
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is less variation across hospitals in the random effects specific to readmission compared to 
the random effects for mortality (either prior to or post-readmission); compare the pos¬ 
terior summaries for SD(V)i) to those of SD(V} 2 ) and SD(V) 2 ). Furthermore, while there 
is no evidence of correlation between hospital-specific random effects for readmission and 
corresponding random effects for mortality, there is some evidence of a positive correlation 
between hospital-specific random effects for mortality pre- and post-readmission, although 
the 95% CIs each cover 0. 

Table 9: Posterior summaries (medians (PM) and 95% credible intervals (Cl)) for standard devia¬ 
tions (SD) of the underlying population distributions for the patient-specific frailties and hospital- 
specific random effects. Estimates of population correlation components, between hospital-specific 
random effects, are also provided. 


Weibull-MVN Weibull-DPM PEM-MVN PEM-DPM 

PM (95% Cl) PM (95% Cl) PM (95% Cl) PM (95% Cl) 


Patient-specific 

frailties 

Ve 

Hospital-specific 
random effects 
SD(Vji) 
SD(V} 2 ) 
SD(H j3 ) 
corr(V)i, V j2 ) 
corr(Vji, V j3 ) 

corr (V j2 , V j3 ) 


1.03 

(0.94, : 

1 .12) 

0.26 

( 0.20, 

0.34) 

0.37 

( 0.28, 

0.47) 

0.37 

( 0.27, 

0.50) 

-0.04 

(-0.40, 

0.33) 

0.06 

(-0.32, 

0.42) 

0.39 

(-0.02, 

0.67) 


1.03 

(0.95, : 

1 .12) 

0.27 

( 0.21, 

0.35) 

0.37 

( 0.28, 

0.47) 

0.37 

( 0.27, 

0.50) 

-0.04 

(-0.40, 

0.33) 

0.06 

(-0.32, 

0.43) 

0.37 

(-0.03, 

0.67) 


0.61 

(0.50, ( 

3.71) 

0.25 

( 0.19, 

0.32) 

0.32 

( 0.25, 

0.41) 

0.33 

( 0.25, 

0.44) 

-0.12 

(-0.44, 

0.23) 

0.03 

(-0.32, 

0.38) 

0.28 

(-0.12, 

0.59) 


0.61 

(0.49, ( 

3.71) 

0.25 

( 0.20, 

0.32) 

0.32 

( 0.25, 

0.42) 

0.33 

( 0.25, 

0.45) 

-0.12 

(-0.45, 

0.24) 

0.03 

(-0.33, 

0.39) 

0.29 

(-0.11, 

0.59) 


6.3.5 Hospital-specific random effects 

As noted, a key advantage of embedding the analysis of cluster-correlated semi-competing 
risks data in the Bayesian framework is the relatively straightforward nature of obtaining 
posterior summaries for the hospital-specific random effects themselves. Figure 0] provides 
posterior medians and 95% CIs for 1 %, j = 1 , ..., 112 , based on the four models in 
which a semi-Markov specification is adopted for /io 3 (•) - Note, across the four panels, the 
ordering of the hospitals is based on the magnitude of the posterior median under the 
Weibull-MVN model. Comparing the panels we see that posterior uncertainty for the V\j 
is generally greater under models that adopt a DPM for the hospital-specific V compared 
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to those that adopt a MVN specification. This may not be surprising given the additional 
complexity of the DPM specification, although we do find that the more ‘complex’ PEM 
specification for the baseline hazard functions yields lower posterior uncertainty than the 
Weibull specification. 

6.3.6 Hospital-specific ranks 

In addition to examining the absolute values of the hospital-specific V\ j, we also considered 
their rank ordering. Figure [5] compares the ranks of the J —112 hospitals according to the 
posterior median of V\ 3 under the PEM-MVN model with a semi-Markov specification for 
h 03 (-) t° the corresponding ranks based on four other models: (a) LN-GLMM; (b) Weibull- 
MVN with a semi-Markov specification for /io 3 (•)5 ( c ) PEM-DPM with a semi-Markov 
specification for h os(-); and, (d) PEM-MVN model with a Markov specification for /io 3 ( - )- 
In each panel, the grey horizontal and vertical lines mark the ‘top 10’ hospitals (i.e. ranks 
1-10) and ‘bottom 10’ hospitals (i.e. ranks 103-112). 

From panel (a) we see that the correspondence between the ranks under a semi-Markov 
PEM-MVN model and those under a LN-GLMM it is far from exact. Crucially, from 
the lower-left portion of the panel, three hospitals that would have been ranked in the 
top 10 under the semi-Markov PEM-MVN model are ranked outside the top 10 under the 
LN-GLMM (specifically, those marked with a ^). Correspondingly there are three hospitals 
(marked with a ▲) who are indicated as being in the top 10 under the LN-GLMM while the 
semi-competing risks analysis under the semi-Markov PEM-MVN would have ranked them 
outside the 10 top. Furthermore, from the top-right portion of the panel, three hospitals 
ranked in the bottom 10 under the semi-Markov PEM-MVN model are ranked above the 
bottom 10 under the LN-GLMM model (i.e. those marked with a A). 

From panels (b)-(d) we find that there is greater correspondence in the ranks of the 112 
hospitals across the models within the proposed hierarchical framework. Comparing the 
ranks under the semi-Markov PEM-MVN specification to the semi-Markov Weibull-MVN 
specification in panel (b) we see that twos hospital that would have been ranked in the top 
10 is now outside the top 10; there is also one hospital that is ranked in the bottom 10 
under the semi-Markov PEM-MVN specification but outside the bottom 10 when a more 
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Figure 4: Posterior summaries (median and 95% credible interval) for the hospital-specific random 
effects for readmission, V±j, under four models for which a semi-Markov specification for hos(-) 
was adopted. In each panel, the hospitals are ordered according to the posterior median under 
the Weibull-MVN model. 
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(b) 



(c) (d) 



Figure 5: Comparison of ranks of J=112 hospitals in the Medicare data on the basis of the 
posterior median for hospital-specific random effects for readmission, V\ y Panels compare the 
ranks on the basis of one of four models to a referent set of ranks based on a PEM-MVN model 
with a semi-Markov specification for hos(-) (see Section f6. 3. 61 for details). Hospitals marked with 
a % suffer under the given model compared to the referent under the referent semi-Markov PEM- 
MVN model; hospitals marked with a A benefit. 
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restrictive Weibull model is used for the baseline hazard functions. Panels (c) and (d) 
are qualitatively similar in that the same two hospitals switch ranks at the lower end and 
the same two switch at the upper end; more generally, consistent with the conclusions we 
draw from Table [7] there is very close correspondence in the ranks across the three models 
represented in these two panels. 

7 Sensitivity Analyses 

As outlined in Section 13.51 and Table [2] the proposed Bayesian framework requires the spec¬ 
ification of a number of hyperparameters. In practice comprehensive sensitivity analyses 
should be conducted to examine the extent to which conclusions are robust with respect 
to this specification, especially across key targets of estimation and inference. Here we 
focus our attention on the choice of hyperparameters for the prior distribution of Ey, 
the variance-covariance matrix of underlying population distributions for hospital-specific 
random effects, and their influence on estimation/inference for the random effects as well 
as r, the precision parameter in the DPM specification of the baseline hazard function. 
Towards this, we conducted sensitivity analyses based on the semi-Markov PEM-MVN 
and PEM-DPM models, specifying a range of values for p v ) and (T 0 , po) such that 
\f/ v ='& Q ='&*(p*-4:)I 3 and p v —po=p*, where, 0*=O.O1, 0.1,1,10 and p*= 5,10, 50,100. Note, 
these specifications correspond to a prior distribution for Ey with a mean of 0*/ 3 and a 
variance of diagonal elements of 20* 2 /(p* — 6). 

Table flOl presents the results. First, we focus Case I-IV, where 0*=1, p* = 100, 50,10, 5 
which correspond to prior distributions of Ey having a mean of J 3 and a standard deviation 
of diagonal elements of 0.15, 0.21, 0.71, 3.16; for Case IV we note that the (induced) prior 
standard deviation was calculated from 100,000 random draws from the prior. From the 
results we see that when the prior distribution is centered around the identity matrix the 
posterior assigns mass to smaller values of SD(V)i) as one increases the prior variance 
(dictated by decreasing p*). This is likely due to the discrepancy between the actual 
variation on the cluster-specific random effects for /qQ and the choice the identity matrix, 
J3, as the prior mean for Ey (since 0*=1) together with the strength attributed to that 
choice (i.e. the prior variance for Ey dictated by p*). When a prior mean of J 3 is chosen for 
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Ey together with a high value of p* the overall prior overcomes the information in the data 
such that the posterior for SD(V)i.) is pushed ‘closer’ to 1 . 0 . As p* decreases, however, and 
less prior mass is given to Ey = I3, the likelihood is able to overcome the less informative 
prior so that the posterior is able to move away from the prior mean. Interestingly, based on 
both the DIC and LPML measures, we find that the overall fit of the data across Cases I-IV 
improves as the prior variance increases. We therefore interpret these results collectively as 
indicating that the variation across the (true underlying) Vj\ is meaningful but relatively 
small. Turning to Cases V and VI, we note that the induced prior distributions of 'Ey are 
centered around relatively small values, specifically O.OI /3 and O.I/ 3 , with induced prior 
standard deviations of diagonal elements of 0.07 and 0.22, respectively. From the DIC and 
LPML values, these specifications of (T*, p*) further improve the overall fit of the model 
from Case IV, we the posterior summaries for SD^x) again indicating that the variation 
in the Vj\ is relatively small. 

While there are clear differences in the posterior summaries for SD(V)i) across Cases 
I-VI, within each case we see that there is little difference in the corresponding summaries 
between the MVN and DPM specifications; that is the conclusions one draws regarding 
the variation of the true underlying Vj\ are robust to this choice. However we do find 
that there are substantial differences in the average posterior standard deviations for the 
J cluster-specific V, 1 . In Case II, for example, ay j , is 0.34 under the MVN specification 
and 0.70 under the DPM specification. Generally, this ordering is consistent across the six 
cases, as well as with the results presented in Figure |4j When combined with the posterior 
summaries for SD(V)i), the results suggest that for our application the trade-off of using the 
more flexible DPM specification is somewhat detrimental to the analyses; use of the DPM 
specification rather than the MVN does not serve to change our conclusions regarding the 
variation across the true Vj 1 but has, rather, served to increase the posterior uncertainty 
regarding any given specific V)i. 

Finally, the last column of Table fTOl presents the posterior median for r, the precision 
parameter in the DPM specification. If one interprets the DPM specification as a mixture 
of MVN distributions (see Supplementary Materials Section B), r dictates, in part, the 
number of mixture components and, hence, the complexity of the overall specification. 
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Table 10: Sensitivity analyses for prior specification for hospital-specific random effects. We obtain 
the fits from PEM-MVN and in PEM-DPM conditioning on a range of values for (T„, p v ) and 
(To, po): T„=To=T*(p* — 4)1 and p v =po=p*, where, ^>*=0.01,0.1,1,10 and p*=5,10,50,100. 
For each model DIC, LPML, and posterior medians (PM) and 95% credible intervals (Cl) for 
standard deviations (SD) of the underlying population distributions for hospital-specific random 
effects are reported, as well as the average posterior standard deviation of hospital-specific random 
effects estimates on readmission (ay jl ). Also presented are the PM for t under the PEM-DPM 
model. 


Case 

(T*, p*) 

Model 

DIC 

LPML 

SD^O 

PM (95% Cl) 

cr Di 

r 

PM 

1 

( 1 , 100 ) 

MVN 

45784.6 

-22905.0 

0.77 

(0.69, 

0.85) 

0.36 


DPM 

45779.3 

-22902.4 

0.77 

(0.70, 

0.85) 

0.70 

0.14 

11 

(1, 50) 

MVN 

45758.2 

-22890.1 

0.66 

(0.59, 

0.75) 

0.34 


DPM 

45754.4 

-22888.0 

0.66 

(0.59, 

0.75) 

0.72 

0.14 

III 

( 1 , 10 ) 

MVN 

45642.5 

-22828.0 

0.39 

(0.33, 

0.47) 

0.27 


DPM 

45644.0 

-22828.7 

0.39 

(0.33, 

0.47) 

0.47 

0.14 

IV 

(1.5) 

MVN 

45574.1 

-22790.9 

0.25 

(0.19, 

0.32) 

0.20 


DPM 

45569.0 

-22789.3 

0.25 

( 0 . 20 , 

0.32) 

0.32 

0.15 

V 

(0.01, 5) 

MVN 

45549.4 

-22777.8 

0.08 

(0.04, 

0.16) 

0.08 


DPM 

45550.6 

-22778.8 

0.09 

(0.04, 

0.18) 

0.12 

0.31 

VI 

(0.1, 5) 

MVN 

45545.7 

-22776.7 

0.14 

(0.09, 

0 . 21 ) 

0.13 


DPM 

45540.6 

-22774.0 

0.15 

( 0 . 10 , 

0.24) 

0.18 

0.23 


From the results, however, we see that the posterior median of r, which takes on values in 
( 0 , oo), tends towards quite small values and is generally robust to the specification of (T*, 
p*). To further investigate the role of r in our analyses, we conducted a series of additional 
analyses where r was fixed at values ranging from 0.1 to 100 (i.e. we did not adopt a gamma 
hyperprior as described in Section 3.5.3). Although details are not reported here, we found 
that results of our analyses to be very robust to the specific value of r, again indicating 
few gains associated with use of the more flexible DPM specification for our application. 
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8 Discussion 


In this paper, we propose a comprehensive, unified Bayesian framework for the analysis 
of cluster-correlated semi-competing risks data. The framework is flexible in that it lets 
researchers take advantage of the numerous benefits afforded by the Bayesian paradigm 
including the natural incorporation of prior information and the straightforward quantifi¬ 
cation of uncertainty for all parameters including hospital-specific random effects. The 
framework is also flexible in that it gives researchers choice in adopting parametric and/or 
semi-parametric specifications for various model components, a key consideration in prac¬ 
tice when small sample size may require pragmatism during the analysis. To facilitate 
model choice, we have also developed DIC and LMPL measures for model comparison 
within the proposed framework. Finally, computationally efficient algorithms have been 
developed and implemented, and are readily-available in a freely-available R package. 

The work in this paper was motivated by an on-going collaboration investigating vari¬ 
ation in risk of readmission following a diagnosis of pancreatic cancer. Towards this, we 
applied the framework to a sample of 5,298 Medicare enrollees diagnosed with pancreatic 
cancer at one of 112 hospitals between 2000-2009. The results from our analysis indi¬ 
cate a number of important determinants of risk of readmission including gender, age, 
co-morbidity status (as measured by the Charlson/Deyo score), whether or not they under 
went a procedure during the index hospitalization, the length of stay of the index hospital¬ 
ization and the location to which the patients was eventually discharged. The analyses also 
revealed that there is substantially less between-hospital variation in risk of readmission 
than the risk of death (either prior to or post-readmission), after accounting for patient 
case-mix. To our knowledge these are the first reported results of this kind in the literature 
and we are currently expanding our analyses to consider patients across the entire U.S. 

More generally, in the clinical and health policy literature, the stand ard an alysis _ap- 


proac. 


1997 


l for inves tigat ing risk of readmission is based on a LN-GLMM ( Normand et ah 


Ashetal. 


2 0121 ). In the specific context of our application, compared with results 


based on the proposed framework, such an analysis yields meaningfully different conclu¬ 
sions regarding which patient-level characteristics are associated with risk of readmission, 
the magnitude and statistical significance of those associations and the ranks of hospitals. 
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Given the relative robustness across models within the proposed framework, the fact that a 
LN-GLMM yields different conclusions is likely related to the fact that death is completely 
ignored as a competing risk. As a concrete example consider the hospital in Figure |5](a) that 
is ranked 8th under the semi-Markov PEM-MVN model and 26th under the LN-GLMM. 
Closer inspection of the raw data reveals that very few patients diagnosed at this hospital 
died within the 90-day window we consider. At other hospitals, the force of mortality is 
stronger and patients die at higher rates within the 90-day window; that these patients die 
is overlooked by the LN-GLMM which assumes that they remain ‘at risk’ to experience a 
readmission event. Hence their estimated readmission rates are too small in the LN-GLMM 
(since the denominator is erroneously inflated). Unfortunately the hospital ranked 8th un¬ 
der the semi-Markov PEM-MVN model suffers from their low mortality rate in the sense 
that they do not benefit from erroneous inflation of the readmission rate denominator, as 
other hospitals do. Hence the change in rank. 

As indicated, results across models within the proposed framework were relatively ro¬ 
bust in our main application. We did find, however, that models which adopted the flexible 
PEM specification for the baseline hazard functions had substantially better fit to the data 
than models that adopted a Weibull specification. While models based on a semi-Markov 
specification for death following readmission generally had better fit to the data than models 
based on a Markov specification we note that this choice does not affect the interpretation 
of the model for readmission (i.e. model (|4]) ) the investigation of which was our primary 
scientific goal. With this in mind, we have not reported on the results for the two models 
for death (i.e. models ([5]) and (|6]) ) although they are available in the Supplementary Mate¬ 
rials E. In practice, researchers may be interested in readmission and death jointly in which 
case t he choice of specification for ho 3 (-) will become critical from a scientific perspective 


(.Lee et ah 


20151 1. In our main application, since the data are relatively rich in terms of 


sample size and the event rates, we have taken the PEM-MVN and PEM-DPM models as 
our primary models for comparison of ranks of hospitals and sensitivity analyses. In other 
less-rich data settings, however, analysts may be in a position where structure is needed 
either in the forms of the baseline hazard functions or for the random effects. Finally, 
we note that in our application a MVN specification for the population distribution of 
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the vector of hospital-specific random effects, V 3l appeared to be adequate. That is, the 
so-called Bayesian non-parametric DPM specification did not yield any additional insight 
into our understanding of variation in risk of readmission nor did it change meaningfully 
the ranking of hospitals. In other applications this may, of course, not be the case and the 
proposed framework gives researchers important choice in this regard. 

In Section 5, we show that incorrect assumption of the underlying distribution for 
cluster-specific random effects or baseline hazard functions result in lower efficiency of 
the incorrect parametric estimators. In addition, the computational efficiency of proposed 
models with non-parametric specification of parameters heavily depends on underlying 
distributions of model parameters. For PEM models, if the underlying hazard function has 
an intricate shape, the model estimates a posterior distribution a g to be centered around 
a larger value, resulting in expensive computation due to more parameters (Ag^’s) to be 
estimated. For DPM models, if data suggest a larger value of r, the model will introduce 
more latent classes in the mixture, implying more parameters to be estimated. 

Our analysis focuses on readmission 90 days post-discharge. However, we note that the 
computational performance of our proposed approach would not be challenged in the cases 
when an administrative censoring is not imposed. In particular, the proposed PEM model 
is flexible in that it allows the time scale for each of three hazard functions to be different 


for each transition. Following iMcKeague and Tighiouartl (2000) and 


Haneuse et al. 


(2008), 


we suggest the last observed event time points be the upper bound in general problems 
where an administrative censoring is not imposed. In our application, however, since most 
of patients diagnosed with pancreatic cancer die within 1-year period, we would expect 
the estimates of baseline hazard functions have a relatively greater uncertainty in the late 
periods if the administrative censoring is not considered. In addition, in the context of our 
study, patients can experience multiple readmission events prior to death. The literature 
on recurrent event semi-competing risks would likely be useful for this setting and thus 
the development of methods that can accommodate recurrent non-terminal events in the 
cluster-correlated data setting is a promising area for future development. 

The proposed hierarchical models assume constant hazard ratios over time conditional 
on the cluster-specific and patient-specific random effects. Since the primary interest of 
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our analysis is the study of readmission event within ‘short time frame’ (30 days or 90 
days after discharge), the proportionality of hazards is quite reasonable assumption in our 
application. In the literature of multi-state models, more rigorous diagnosis can be most 


naturally done by considering a more flexible multi- state rno c 
or inclusion of non-proportional covariate effects Ihlougaard . 


cl such as a stratified model 


20001 ). Expanding the scope 


of the proposed models to include deviation from proportional hazards as well as time- 
varying covariates is our future work. In this paper, we considered a gamma distribution 
for the within-patient frailty because of its computational tractability. When the frailty 
distribution is misspecihed, the resulting estimator is not guaranteed to be consistent, with 
the extent of asymptotic bias depending on the d iscrepancy between the assumed and true 


frailties distributions. However, 


Hsu et al. 


(120071 ) studied the effect of mis-specihcation of 


frailty distribution on the marginal regression estimates and hazard functions when gamma 
distribution is assumed. Their results show that the biases are generally low, even when 
the true frailty distribution is substantially different from the assumed gamma distribution. 
Therefore, if the regression parameters and hazard function are of primary interest, the 
gamma frailty model can be a reasonable choice in practice. 

Finally, we conclude by emphasizing that the proposed framework significantly improves 
and expands the set of statistical tools researches have to study quality of end-of-lifc care. 
While our focus has been on pancreatic cancer, the proposed framework is broadly appli¬ 
cable to all ‘advanced’ health conditions for which current treatment options are limited 
and the force of mortality is strong. Such studies will be of paramount importance in 
the near-future because many of these conditions, including other cancers as well as neu- 
rodegenerative conditions such as Alzheimers’ disease, directly affect large segments of an 
increasingly aging population. In addition, although it has not been in the focus of this pa¬ 
per, the proposed framework will also be critical in helping policy-makers understand and 
ultimately control the increasing costs of health care delivery in the U.S. In particular, the 
proposed framework provides CMS appropriate statistical tools with which to expand the 
scope of the Hospital Inpatient Quality Reporting Program and the Readmission Reduction 
Program to include to conditions with strong forces of mortality. 


SUPPLEMENTARY MATERIALS 
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Title: In online Supplementary Materials, we provide a detail description of Metropolis- 
Hastings-Green algorithm to fit our proposed models. Additional details regarding 
the Medicare data and results from the application are also provided, (pdf hie) 

R-package ‘SemiCompRisks’: R-package SeraiCorapRisks contains codes to implement pro¬ 
posed Bayesian framework described in the article. The package is currently available 
in CRAN. 
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Introduction 


This document supplements the main paper titled “Hierarchical models for clustered 
semi-competing risks data with application to pancreatic cancer”. In Section A, we provide 
technical details regarding the MVN-ICAR specification for baseline hazard functions in 
PEM models. In Section B, we provide a detailed description of the Metropolis-Hastings- 
Green algorithm to implement our proposed Bayesian framework (Weibull-MVN, Weibull- 
DPM, PEM-MVN, PEM-DPM). In Section C, we examine the potential use of methods 
proposed in Gorhne and Hsu (2011) and Liquet et al. (2012) in the context of the mo¬ 
tivating application. In Section D, we provide results from simulation studies that were 
not presented in the main paper. Finally, Section E supplements our main paper with 
some additional results from analyses of Medicare data from New England and a visual 
assessment of convergence of the proposed MCMC schemes using potential scale reduction 
factor. 

In order to distinguish the two documents, alpha-numeric labels are used for sections, 
tables, figures, and equations in this document while numeric labels are used in the main 
paper. 
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A MVN-ICAR Specification for 


In MVN-ICAR, the specification of a prior for the components of X g is considered as a one- 
dimensional spatial problem. The dependence between neighboring intervals are modeled 
via a Gaussian intrinsic conditional autoregression (ICAR) (Besag and Kooperberg, 1995). 
Let A g ^ denote the vector given by X g with the k th element removed. The full conditional 
prior for X gyk is then taken to be the following normal distribution: 


Ag, fc |A^ k) ~ Normal^, a 2 gjk ), (1) 

where the conditional mean, u g)k = ^\ g + )Tb/fc W k] (— p Ag ), is a marginal mean plus a 
weighted sum of the deviations of the remaining intervals. Let A 9 k = s g)k — s g ,k- 1 denote 
the length of the I 9yk interval. We determine the weights for the intervals adjacent to the 
k th intervals based on these lengths as follows: 


w a _ _ c v ( A fc-i + Ap 
Hk-D a ? ,+2Al + Ai 




W 9 = _ C A g (A 3 k + A 9 k+1 ) 

k(k+ 1) A 9 i i A 9 


( 2 ) 


where the constant c\ g G [0,1] dictates the extent to which A 9i fc is influenced by adjacent 
intervals (Haneuse et ah, 2008). The remaining weights corresponding to intervals which 
are not directly adjacent to the k th interval are set to zero. The conditional variance cA k 
in (1) is given by cr\ g Q 9 k . The o\ g is an overall measure of variation across the elements of 
X g and the diagonal matrix Q 9 k is given by 

At, + 2A 9 k + A 9 k+1 ' 1 J 


Given (1), (2), and (3), we can see that X g jointly follows a (K g +1)- dimensional multivariate 
normal (MVN) distribution: 


MVN A ' 9+1 (/i Ag l, <S Ag ), (4) 

where /i Ag is the overall (marginal) mean, a Ag the overall variability in elements of X g . The 
S Ag is given by (/ — W 9 )~ 1 Q 9 , where a ( K g + 1) x (K g + 1) matrix W 9 k j)—W kj and a 
(Kg + 1) x (K g + 1) diagonal matrix Q g {k)k )=Q 9 k - 
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B Metropolis-Hastings-Green Algorithm 


B.l Weibull models 

Let <3 >vk = {a Wj i,a Wt 2 ,a Wt i, K Wt 2 , k w , 3 , fli, 02, /33,7,V} be a set of parameters in the 
likelihood function of Weibull models. The observed data likelihood Lw(P\$w) is given 
by 

J n i 


1 — |-t —p / aw 1 — \ \ h*i(1 h*2) / „ a w 1 — 1 ctiu 3 —1 

( 7?i®tu,l ^wXUjil Vjil J ( 11jil(^w,3^w,3yji2 Vji 3 


hi 1 h?2 


J = 1 i=l 


a ^\ h?-2 (^ hil ) 

X \ r )jiOi w 2^w,2'!jji2 Vji2J exp { T'W (j/jilj Uji2) / j 


( 5 ) 


where = exp (a?!^ + V ig ) and 

r w{tjil, tji 2 J 

lj% {^w,it a ji{ lr ljii + K w rf]i\r]ji 2 + (k w , 3 t^' 3 ~ K w, 2 t°ji\ 3 ) Vjtt} , for Markov model 
7 ji [ K w,it < ja 1 Vjii + ^ w , 2 t a ji{ 2 r]ji 2 + {n w ,3{tji2 - tjn) aw ’ 3 } r) ji3 ] , for semi-Markov model 
For Weibull models, we use a random scan Gibbs sampling scheme, randomly selecting 
and updating a (vector of) model parameter at each iteration. 


B.1.1 Updating (3 g 

Let denote a set of parameters <f> with (3 removed. The full conditional posterior 

distribution of /3i can be obtained by 

Tr^l^UEy) cx L W {D\$ W ). 

J rij 

oc IIII exp _ 7ji«t U ,i2/5i ,1 e a! ^ 1 ^ 1+Vj ' 1 

i=i i=i 

Analogously, the full conditionals of /3 2 and /3 3 are given by 

J rij 

7T(/3 2 |$ w (/32) , 0, £y) oc JJ n exp |^ 2 (1 - 5jn)xJ i2 f32 - ^jin W) 2 y^i’ 2 e x ^ 2+Vj2 1, 

J=1 i =1 
J 

7r (/ 3 3|^w^ 3) , S v) OC JJ II exp - 7ji«u;,3 ( 2 /S 2 3 - | 

j=l«=1 
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Since the full conditionals do not have standard forms, we use Metropolis Hastings (MH) 
algorithm to update each element of (3 g , /3 9) i, ... ,/3 g , pi . In our algorithm, the conventional 
random walk MH is improved in convergence speed by taking some meaningful function of 
the current value /3^ k ' 1 for the mean and variance of Normal proposal density. Specifically, 
let Di(f3 g: k) and D 2 (/3 9t k) denote the first and second gradients of log-full conditional of 
f3 g with respect to f3 g ^, then a proposal j3* is drawn from a Normal proposal density 
that is centered at /i(/3^' fc '^) = fj^' k l> — Di(/3g k l ^) / D 2 (j3^ k ^), the updated value from the 
Newton-Raphson algorithm, with a variance cr 2 (/3^ k 1 ' > ) = — 2A 2 / D 2 (/3g k 1 ' 1 ), based on the 
inverse Fisher information evaluated with /3 gk (Roberts et ah, 2001; Gelrnan et ah, 2013). 
Therefore, the acceptance probability for /3 g> k is given by 

6 ,EyJNormal 

7 r(/ 3 < , - 1 ) | 4 - ( ' , »>, 0, £„)NormaI (’ 

where {3^ 1 is a sample of (3 g at current iteration and (3 * is the (3 g with fc-tli element 
replaced by j3*. 


B.1.2 Updating a w>g 

The full conditional posterior distribution of a w> i is given by 

Tr{a w ^ aw ' l \d^ v ) 
oc L W (D\$ W ) x Tr(a wA ) 


J rij 


oc i 1 l e 6q ’ iQu) ’ 1 Yl 11 ex P (iji^iy^Vjii) ■ 

j =1 i =1 


Analogously, the full conditionals of a W)2 and a w ^ are given by 

J rij 

/ I af. (^ia,2) /l Y" 1 \ f^a,2 1 —brr lOLin 9 T T T T I &w,2 \ ^ji2 (1 h'il) / CX W 2 \ 

oc a w ’ 2 e ’ ’ [[ ( a w,2y ji2 ) exp {-'y j iK w , 2 y jil 'Vju) , 

j =i i=i 
J fij 

/ I af. (^ia,3) /I Y" 1 \ 3 1 —brv 30 i ll , 9 "1 T "1 T ( &w,3 \ 

Ot a w,3 e ’ ’ \ a w,3yji 2 ) 

j =1 i =1 

x exp |- 7 jiK W) z(y^ - Vj^Vjis] ■ 

In MH algorithm to update a W:9 , we generate a proposal a* from a Gamma distribution 
with Gamma ((a™,/' 1 ) 2 //^, ) which corresponds to a distribution with a mean of 
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a>w,g and a variance of fc 0 . The value of ko is specified such that the MH step for a Wt9 
achieves an acceptance rate of 25% ~ 30%. Finally the acceptance probability to update 
a w , g can be written as 

K{a*\^ aw '«\e,Y, v )g (a& 1} |« s ) 2 /fc 0 ,< g Ao) 

\^ aw ’ 3 \0, Zv)G (alj^-y/ko^^/ko) ’ 

B.1.3 Updating k W}9 

The full conditional posterior distribution of K Wt9 can be obtained by 

7r(K W! g\^ Kw,B \d, Sy) oc L W (D\$ W ) x n(n Wt g). 

We see that the full conditionals of k, W)9 are gamma distributions and the samples can be 
drawn from following distributions: 


k w , 2 \^ Kw ’ 2) ,6,Z v 
Kw,3\^w Kw ' 3 \Q-, Sy 

B.1.4 Updating 7 ^ 

The full conditional posterior distribution of 7 ^ is given by 

oc L W {D\$ W ) x n(^ji\6) 

oc 7 ji exp [-r w (y jil: y ji2 ) - 0 TnJ ■ 

Therefore, we sample 7 ^ from 

Gamma (Sjn + S ji2 + 6 > _1 , 1 , y ji2 \ 77 = 1) + 6 1 " 1 ) . 


( J n j J n j \ 

EE 5jil + 0/c,l5 EE + 6*,i I, 

j=l i=l j=l j=l / 

( J % J nj \ 

EE ^2(1 — $jil) + 0/t,2, EE IjiVjil Vji2 E b K) 2 j , 

J=1 i=l j=l i=l / 

( J rij J rij 

EE bjnSji 2 + a Kj 3, EE in (vji7 3 - y% r 3 ) vjis + &*,a 

j=l «=1 j=l i=l 
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B.1.5 Updating 9 


Let ^ — 1/6 denote the precision parameter of frailty distribution. The full conditional 
posterior distribution of £ is given by 


J n 3 


7r(£|$iv,£y) oc tt(£) 

j =1 i =1 


OC 


Cnt+b g -l „-£(E/=i Eiii 7 ji+o-e) 

' imr - 


J rij 


1111 ; 


e-1 


3= 1i=1 


We revise the traditional random walk MH algorithm for updating £ as done in Section 
B.1.1 for /3 g . Let //$(£) = £ - min{0, L>i i? (£)/U 2 , ? (£)} and cr|(£) = -c 0 /D 2> ^), where 
U!^(£) and U 2 ^(£) are the hrst and second gradients of log7r(£|<f)^^, Ey) with respect to 
£. A proposal £* is generated from the following Gamma distribution 


Gamma , J K ( ‘“ ,, )/<’?K““ 1> )) • 


The value of Co > 0 is specihed such that the algorithm achieve the desired acceptance rate. 
The acceptance probability to update £ is then given by 


7T(£*|<fy, Ey) Gamma (£1^(D 2 /^ 2 (D, MD/^f (D) 
7r(£h _1) |$iy, Ey) Gamma (£*|/x^(£ (t-1) ) 2 /cr|(£ (t - 1 )), /x(£ (t_1) )/cr|(£G” 1 ))) 


B.1.6 Updating V) for Weibull-MVN model 

The full conditional posterior distribution of Vj\ can be obtained by 

'K{V jl \^ Vjl \e^v) oc L W (D\Q W ) x 7r("V£|Ey). 

« ex P jxj ( v n s 3ii - 

Analogously, the full conditionals of Vj 2 and V) 3 can be written as 
n(V j2 \^ 2) ,9, E v ) (x exp j (V j2 5 ji2 (l - S ja ) - - ^ T S ^Vj \ , 


i =1 

Tlj 


<V j3 \^ Via) , e : Ey) oc exp <j ^ (V j3 5 ja s ji2 - ~ ^V/E^V- \ . 


i=l 


As done in Section B.1.1, in a MH step for updating Vj g , we sample a proposal V* from 
a Normal distribution that is centered at Hv(Vjl ~= Vjg~^ — Diy{Vj y ^ 1 ' 1 )/D 2 y{Vj t g _1) ) 
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and has a variance of = —2A 2 /D2y(V < ' t ~ 1 ' 1 ), where D 1 y(Vj g ) and D 2 y{Vj g ) are 

the first and the second gradients of log7r(V 7S | < f ) w } V " J9 \ 0, Ey) with respect to Vj g . Finally, 
the acceptance probability is given by 

-K(V*\^ Vj3 \ 6 , Ey) Normal (y^% v (V*), a 2 v (V*)^ 

n(v%- 1) \^ Vja \o,z v )rnrm a \ (V*|M ^)’’ 

B.1.7 Updating Ey for Weibull-MVN model 
The full conditional posterior distribution of Ey can be written as 

j 

7r(Ey|<f>ty, 0) OC 7r(Ey ) 7r( Vj | Ey ) 

3 = 1 

« IWI-^expj-i^^ + vl? 

Therefore, we update Ey from the following inverse-Wishart distribution: 

Ey|$vv, ^ ~ inverse-Wishart E v > v j j + J +p< 

\3 =1 

B.1.8 Updating Vj and Ey for Weibull-DPM model 

Towards developing this model, suppose that, instead of arising from a single distribution, 
the Vj are draws from a finite mixture of M multivariate Normal distributions, each with 
their own mean vector and variance-covariance matrix, (/x m , E m ) for rn — 1,..., M. Let 
rrij e {1 denote the specific component or class to which the j th hospital be¬ 

longs. Since the class-specific (/x m , E m ) are not known they are taken to be draws from 
some distribution, G q. Furthermore, since the ‘true’ class memberships are not known, 
we denote the probability that the j th hospital belongs to any given class by the vector 
p = (pi, ■ ■ ■ ,Pm) whose components add up to 1.0. In the absence of prior knowledge 
regarding the distribution of class memberships for the J hospitals across the M classes, a 
natural prior for p is the conjugate symmetric Dirichlet(r/M,... ,r/M) distribution; the 
hyperparameter, r (Walker and Mallick, 1997). Jointly, this finite mixture distribution can 






be summarized by: 


Vjlrrij ~ MVN ( p m ., E m .), 

Om, S TO ) ~ G 0 , for m = 1,..., M, 
mj\p Discrete(mj| pi,... ,Pm), 

p ~ Dirichlet (r/M,..., t/M), (7) 

Finally, letting M —y oo the resulting specification is referred to as a Dirichlet process 
mixture of multivariate Normal distributions (DPM-MVN) (Ferguson, 1973; Bush and 
MacEachern, 1996). When M —y oo, we cannot explicitly represent the infinite num¬ 
ber of (p m , £ m ). Instead, following Neal (2000), we represent and implement the MCMC 
sampling for only those (p m , £ m ) that are currently associated with some observations at 
each iteration. In this subsection, we provide a step-by-step detailed description of the MH 
algorithm to update Vj in Weibull-DPM model. 

First, we update a class membership rn 3 based on mj\rri(-j), Vj, j — 1, • ■ ■ , J. Let m_(j) 
denote a set of all class memberships from clusters except the cluster j. After identifying 
the “n m ” unique classes of m_u\, we compute the following probabilities for each of the 
unique values m. 

P(m j = m\m(- j ),Vj) = J Normal(V)| p m ., E m .)dH_ jiTn (p, £)(8) 

P(rrij ^ m k , \/k ^ j \m^j),Vj) = h j + T j Normal (V)|/z, S)dG 0 (/i, £), (9) 

where H_ 3m is the posterior distribution of (/i, £) based on the prior Go and {Vk : k ^ 
j,rrik = c}. The normalizing constant b makes “n m + 1” probabilities above sum to 1. Let 
A = {j : rrij = m} and Ha be the posterior distribution of (p, a) based on the prior G 0 and 
{Vj : j G A}. It can be shown that the Ha is also Normal-inverse Wishart distribution as 
G 0 is conjugate to multivariate normal distribution: 

1. we draw a sample of a class membership. 

i) For each rrij, identify the n m unique values of m 

ii) For each of the unique values m, compute the following probabilities: 


P{rrij = m Vj) 

u n ~h™ 

f N011118,1 ( Vj | fAjrrij ? ^ rrij ) 


b J-l + rJ 




(10) 

Pinij ^ m k , Vk ^ j |m ( _j), Vj) 

- h T 

f Normal (Vj/Lt, H)dGo{p, £), 

(11) 

1 

5 

1 

+ 


9 






where H_^ rn is the posterior distribution of (//,£) based on the prior Go and {V 4 : 
k 7 ^ j, nik = m}. The normalizing constant b makes n m + 1 probabilities above sum 
to 1.0. Let A = {j : mj = m} and H a be the posterior distribution of (/i, a) based on 
the prior Go and {Vj : j G A}. It can be shown that the Ha is also Normal-inverse 
Wishart distribution as Go is conjugate to multivariate normal distribution: 

H a (h, T,\h a ,Ca, Va,Pa), 


where 


Pa = 


5M0 + W4 


E + W 


Co 


-1 


5 (a — \ -r + |kl| , Pa — Po + \A\, 


M 

Co 


*x = *.+S3 (n - v A ) (Vj -v A y+ -t-tjxi ^ ^' i(12) 

Co ' ' 


jeA 


with V A — m J2k&A V k- Now we define 


Q(Vj, ho, Co, T 0 , Po) 

/ /AC 3 (^j|h, S )dF N iw{p , 2|/io, Co, T 0 ,po) 


IT I P0 

T n “ 


To + VjV T + i-hohci" 


1 + * 


-1 


i/x 0 + VC 




V,- 


PQ + 1 


x- 




X 


rush 

r„,3(f) 


(13) 


It follows that the integrals in (10) and (11) are equal to Q(Vj, pi a, (a, Ta, Pa) and 
Q(Vj, fjL 0 , Co, T 0 , po), respectively. 


iii) Sample m^ ncw ' ) 


based on the probabilities given in ( 10 ) and ( 11 ). 


2. For all m G update (p m , S m ) using the posterior distribution that is 

based on {Vj : j G {k : rrik = m}}. 

3. For j = 1 update Vj using its full conditional using Metropolis-Hastings algo¬ 

rithm. 

4. We treat r as random and assign gamma prior Gamma(a r , b T ) for r. Following Escobar 
and West (1995), we update r by 


i) sampling an c G (0,1) from Beta(r + 1, J), 
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ii) sampling the new r from the mixture of two gamma distributions: 

p c Gamma(a r + n m , b T — log(c)) + (1 — p c )Gamma(a r + n m — 1, b T — log(c)), 
where the weight p c is defined such that p c /(l — p c ) = (a T + n m — l)/{J(b T — log(c))}. 
5. Finally we calculate the total variance-covariance matrix: 

1 J T 

Sy = 'y ^ A) (a 4 m,j A) “1“ } j (14) 

i=i 

where A = 1 VmJJ- 
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B.2 PEM models 


Let $p = {Ai, A 2 , A 3 , /3 i, /3 2 , 03, 7 , V} a set of parameters in the likelihood function of 
PEM models. The observed data likelihood Lp(T>\$p) is given by 


J n i 

nn 

j =1 i= 1 


'K i+l 


lljiVjn exp < ^ Aifc/(si )fe _ 1 < y ja < s ltk ) 


k =1 




X 


X 


rxi+1 ie 3 +i 

l‘ji r 1jil r 1ji3 ex P < ^ Aifc/(si jfc _i < Uji 1 < Si,fc) + ^ A 3 fc/(s 3i fc_i < y ]l2 < s 3 ,fc) 
l fc=l fc=l 

/" P'2 + 1 ^ 1 — 




lfjiVji2 exp < A 2 fc/(s 2 ,fc-i < Vja < s 2 , fe ) 
l fc=i 

x exp {—r'p(yj,;i, yp 2 )} , 


(15) 


where r]j ig = exp ( xj ig /3 g + V jg ) and 


Tp(tjj {, tjj 2 ) 


[ 7* (»i,i Efif + +« Efit' <*‘A?» 

{+, (»;,! Ehd A‘ jt + , jj2 Efid ^.‘A^ 


+ 1*3 ES 1 ^*^), 

+ » 2 i3E£ke A «A-f l ), 


for Markov model 
for semi-Markov model 


= max 


A*f z = 


jo, min(i/ ja , s Sifc ) - s 9 ,fc-i}, 

max jo, min(y ii2 , s 9> ;) - max(|/ ifl , s 9 ,/-i) j>, 
max jo, min(y Jj2 - y ja , s g j) - %j_i)j, 


B.2.1 Reversible jump MCMC algorithm 


for Markov model, 
for semi-Markov model. 


For PEM models, we use a random scan Gibbs sampling scheme, randomly selecting and 
updating a (vector of) model parameter at each iteration. Let BI g and DI g denote a birth 
and a death of a new time split for transition g G (1, 2, 3}. The probabilities for the update 
BI 9 and DI 9 are given by 


71 


K a 

Bln 


7T 


K g 

DI a 


f Poisson (A" g + l|a^„) "i 

p n mm < 1 , ---;- p— > 

l Poisson(A 9 |«p' g ) 1 

r Poisson (K g - l\ctK g ) 1 

p u nun < 1 , ---;-—a- \ 

l Poisson(A 9 |«p' g ) / 


• li aK o 1 

p a mm < 1, ->, 

Ha l ’ K g + ir 

p g min 11 , — ), 

1 aK g 
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where p g is set such that n + 7r^f < C 9 and < 1 for A" 9 = 1,..., K g max . 

K g inax is the preassigned upper limit on the number of time splits for transition g and we 
set 7r^ max = 0. The probabilities of updating other parameters are equally specified from 
remaining probability 1 - + ^Di g )- 


B.2.2 Updating (3 g 

The full conditional posterior distribution of (3\ can be obtained by 

7r(/3i|$p (/3l) ,/x A ,cr^0,Ey) 
oc Lp(D\$ P ) 

J rij / K x + \ 

°C nn-P fijiixjufli - 7 ■ ji e x » lfh+v i 1 eXl,kA )ik I 1 


j= 1 i =1 
\T 


k =1 


where [i\ = (p\ 1 , p\ 2 , p\ 3 ) 1 and er A = (cr Al , o\ 2 , cr| 3 ) T . Analogously, the full conditionals of 
(3 2 and (3% are given by 


J Uj 


K 2 + 1 


oc nn-P ^<2(1 - Sjii)x] i 2 ( 3 2 - ~ji< A-*’' 1 -' •' J ^2 e A2 ' ; A 
3= 1 *=1 ' 

,-(Ps) .. _2 


jil ( ’ 


«=1 


7r(/3 3 |$ P l ,/i A , <r A , 0,1! 


vj 


J n j 


K 3 + l 


°C nn-P - 'y ji e x ^ +Vi3 ^ e X 3 ’ m A$ m J . 

j=l i=l \ m=l / 

As the full conditionals do not have standard forms, we use MH algorithm to update 
each element of (3 g . A detailed description of the adapted random walk MH algorithm is 
provided in Section B.1.1. 


B.2.3 Updating \ g 

The full conditional posterior distribution of Ai is given by 
vr(Ai|$p (Al) ,/x A ,cr2,0,Ey) 
oc L P (U|$p) 7 r(A 1 | / u Al ,cr Al ) 

J nj 

OC nn exp {8jii\i k I(s hk -i < Ujii < si,fc) - 7pAh fc e Al %a} 

j =l i=l 

X 6XP {~2 ^~^ Al ~~ ^ 1 1 ) T S Ai 1 (A i - /i Al l| , 
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where 1 denotes a K g + 1 dimensional vector of l’s. Analogously, the full conditionals of 
A 2 and A 3 are given by 



J n i 


oc nn-P 2(1 - S jil )X 2k I(s 2jk - 1 < y ji2 < s 2 ,k) - 7 jA%k eX2k Vji 2 } 


3= 1 i=l 



^■(A 3 | *;<*•>,9, e v ) 


OC nn exp {^a^i 2 A 3 fe/(s 3 ,fe-i < 


^ Uji 2 A S37) 3/c ^i 3 J' 


i=i i=l 


x exp -^(A 3 - /r A 3 l) T S A 1 (A 3 



— (^3 — hA 3 l) T S A3 1 (A 3 — H\ 3 1) 


Since the full conditionals do not follow known distributions, MH algorithm is used to 
update each element of X g . We follow the adapted random walk MH algorithm described 
in Section B.1.1. 

B.2.4 Updating 7 

The full conditional posterior distribution of 7 ^ is given by 


oc Lp(D\$ P ) x TT(jji\ 9 ) 



exp [-r P (y jil ,y ji2 ) - 9 1 7 J - i ] . 


Therefore, we sample 7 ^ from 


Gamma (S jn + S ji2 + 9 1 , r P (y jil , y ji2 ; 7 ^ = 1) + 9 x ) . 


B.2.5 Updating (y g , a g ) 

Full conditional posterior distributions for y\ g and v g = l/cr A9 , g = 1,2,3 are Normal 
and Gamma distribution, respectively. Therefore, we use Gibbs sampling to update the 
parameters. We obtain the posterior samples of y\ g from 
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because the full conditional is given by 




1 T S^ 1 1 

K eX P ^ 2a 2 I PA„ - 


l^A, 

1 t E7'1 

A 9 


We update v g = 1 /cr A 2 from a Gamma distribution given by 

K g + 1 


Gamma a a , g + 


+ 2 ’ g ^ ^s) ^Aglp-Agl Aj,)^ , 


as the full conditional of v g is 


vr(n 9 |<f>p,/x A , (<r 


2 ^ {<t ' 9 ) ,0,E 


vJ 


« 7r(A 9 |/i A9 ,a^)7r(n g ) 

, Kg +1 1 

°c (n s ) aCT ’ 9+ 2 exp 


&(J,g + 2 (hAg 1 A g ) S Atj (/i As l AgJ f Vg 


B.2.6 Updating 6 


Updating the precision parameter £ = 1/9 in PEM models requires the exactly same step 
as that in Weibull models. Therefore, the readers are referred to Section B.1.5 for the full 
conditional posterior distribution of £ and the MH algorithm. 


B.2.7 Updating Vj for PEM-MVN model 

The full conditional posterior distribution of Vj\ can be obtained by 

'K{V jl \^- {Vjl \^crl9^ v ) oc Lp(D\$p) x 7r(^|E y ). 


K i + l 


OC 


exp { ^ ( V^Sja - 'Yjirjjn ^ e Xl - k A) ik 


i =1 


k =1 


Analogously, the full conditionals of V£ 2 and Vjs can be written as 


’rfel* p (V "\ Ma, o-t 9, Sv) « «pT( U^W 1 - S ja ) - 7ji%a E ^ A i« ) “ }■, 


iC 2 +l 


2 > • 


2=1 

Tin 


1=1 


K 3 +I 


7r(U i3 | < f ) P ( ^ 3) , p A , cr 2 , 0, E y ) oc exp <{ ( V^SjuSja - IjiVjis J2 ) “ 2 V j ' V i 


c , *3 I T/T^-l 

2=1 \ 772=1 

For updating each element of Vj, we use the adapted random walk MH algorithm provided 
in Section B.1.6. 
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B.2.8 Updating Ey for PEM-MVN model 


The full conditional posterior distribution of Ey in PEM-MVN model is the exactly same as 
that in Weibull-MVN model. Readers are referred to Section B.1.7 for the Gibbs sampling 
step for updating Ey. 


B.2.9 Updating V 3 and Ey for PEM-DPM model 

Updating Vj and Ey in PEM-DPM requires the exactly same step as that in Weibull-DPM. 
Therefore, the readers are referred to Section B.1.8 for detailed algorithm to update Vj and 
Ey. Note that in step 3 of the algorithm, the full conditional of Vj needs to be obtained 
based on L P (D |$p) for PEM-DPM. 


B.2.10 Update BI 

We specify log /i 0 (t) = Y^k^i 1 -I(t e I g ^) for the baseline hazard function corresponding 
to transition g with partition (K g , s g ). Updating ( K g , s g ) requires generating a proposal 
partition and then deciding whether or not to accept the proposal. For update BI (a birth 
move), we first select a proposal split time s* uniformly from among the observed event 
times which are not included in the current partition. Suppose s* lies between the (k — l) th 
and k th split times in the current partition. The proposal partition is then defined as 

(0 ■■■) Sg,k— I; S ) ^g,ki ^g,K i+l ^g, max) 

_ / r\ * * * * * _ \ 

= ~ S g, O’ ■■■) S g,k -15 S g,k+ 1; ') S g,K 1 +2 ~ S g, maxj- 

A height of the two new intervals created by the split at time s* also needs to be proposed. 
In order to make the old height be a compromise of the two new heights, the former is 
taken to be the weighted mean of the latter on the log scale: 


( S S g,k-l)^g,k A ( S 9,fc S )^g,k +1 ~ ( S 9,fc s g,k-l)^g,k- 


Defining the multiplicative perturbation exp(A* J+1 )/exp(A* ■) = 
Uniform(0,1), the new heights are given by 



A 


g,k 


$9 ' k _ S * log 

s g,k S g,k- 1 \ U J 


(1 - U)/U, where U ~ 
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and 


1 - U 
U 


A 


* 

9,k +1 


A g^k T 


S Sg,k— 1 
s g,k ~ s g,k -1 



The acceptance probability in the Metropolis-Hastings-Green step can be written as the 
product of the likelihood ratio, prior ratio, proposal ratio, and Jacobian. For g — 1, they 
are given by 


likelihood ratio 
prior ratio 


proposal ratio 


Jacobian 


L P (D\$*p) 

L P (D\$ P y 

Poisson (Ad + l\a Kl ) x MVNk 1+ 2 (A*|/x Ai 1 , o^E^) 
Poisson(AJ| a Kl ) x MVN* 1 + i(Ai|// Ai 1, o^EaJ 
x (2Ad + 3)(2Ai + 2)(s* — si,fc-i)(si,fc — s*) 

s l,max( s l ,k — s l,k-l) 

^DI X (1 /(Kl + 1)) 

Kbi X (1/1 \{y ja : Sji 1 = 1}) x Uniform(f/|0,1) 
pmin(l, (Ad + l)/Q Xl )ft{|/iii : Sji i = 1} _ ttfaii : Jpi = 1} 
pmin(l, 0^/(1 + A'i))(AJ + 1) a Kl 

d\lk/d^i,k d\\ k /dU 
d^i : k+i/d^i,k d\ lk+1 / dU 


u{i-uy 


(16) 


where $p is $p with Ai replaced by A*. 


B.2.11 Update DI 

For update DI (a death or reverse move), we first sample one of the K g split times, s 9t k- 
The proposal for time splits is given by 

(0 Sg t o, •••, Sg k —\^ ^g,k- 1-1, •••> Sg,K g +1 $ g ,max) 

_ / /~V _ =}= 5}C 5(C _ \ 

= — s g, 0; S g,k- 1) •••) S g,K g — s g,maxJ- 


Following Green (1995): 


($g,k Sg,k— l)Ag i fc T (s^fc+l Sg : k) A^fc+l (>Sgr,fc+l ^g,k—l)Xg k , 

e A s .fc+i 1 — t/* 

perturbation : —t-= ———, 

e A S’ fc A* 
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where U ~ Uniform(0,1). Then the acceptance probability can be obtained as the product 
of following four components (for g — 1): 


likelihood ratio 
prior ratio 


proposal ratio 


Jacobian 


L P (D\$*p) 

L P {D\$ P y 

Poisson(ifi - 11 a Kl ) x MVN^A^I/ia,!, a Al S Al ) 
Poisson(iPi|« i ^ 1 ) x MVN A - 1+ i(Ai|/i Al l, cx^Eai) 

S l,max( S l,fc+l — s l,fc-l) 

X_:_ 

(2-ffd T Si,fc— i)( si,fc+i ^i,fc) 

TIBI X (1 /i{yjii ■■ Sjn = 1}) 

Td/ x (1/-^i) 

p min(l, OiKx/Ki)Ki _ a Kl 

pmin(l, K ]_/: 5 ja = 1} Jj {y ja : 8 ja = 1}’ 

d^ik/d\ lk d\\k/d\ l k+1 _ JJ(1 U ) 

dU / dA* fe dU/ d \\ k+1 
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C The potential use of existing methods 


The methods in the main manuscript were developed specifically for on-going collaboration 
examining the risk of readmission following a diagnosis of pancreatic cancer. As indicated 
in the manuscript, the current standard for the analysis of cluster-correlated readmission 
data is a logisitic-Normal generalized linear mixed model. This model ignores death as a 
competing risk, however, and, as such, is inappropriate in for the study of pancreatic cancer 
due to its strong force of mortality. 

Viewing the data arising in the pancreatic cancer as cluster-correlated semi-competing 
risks data , the existing literature does have a number of options that could be considered. 
Here we review two of these options, specifically those proposed in Liquet et al. (2012) and 
Gorfine and Hsu (2011). For the former, we note that the methods have been implemented 
in the frailtypack package for R. 

For convenience, expressions (4)-(6) from the main manuscript that the describe the 
key features of the proposed hierarchical model are repeated here: 

= jji h m {tjn) exp-fX^/?! + V^}, t ja > 0 

Xji2,Vj2) = 7 ji ho 2 (tji 2 ) exp {Xj i2 f3 2 + V j2 }, tj i2 > 0 

hz(tji 2 \tji\, Xj i3 , Vjs) = '-fji h 03 (tj i2 \tjii) exp{Xj i3 f3 3 + V} 3 }, 

tji 2 tji 1) 

C.l Liquet et al. (2012) 

The R package frailtypack provides several classes of frailty models for multivariate 
survival data including shared frailty models, additive frailty models, nested frailty models, 
joint frailty models (Rondeau et ah, 2012). Among these, the shared frailty model and the 
joint frailty model are most relevant the context we consider; additionally, these models 
form the basis for the analyses presented in Liquet et al. (2012). Here we provide a summary 
of these two classes using the notation developed in the manuscript, as well as an overview 
of their drawbacks in regard to the analysis of cluster-correlated semi-competing risks data. 
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C.1.1 The shared frailty model 

In the shared frailty model, the hazard function for the subject i in the cluster j conditional 
on the cluster-specific shared frailty term rjj = (rjji, rjj 2 , rjj 3 ) is given by 

= rjjih Q , (t^ j) exp{ Xj t ^i}, t ja > 0 

h 2 {t ]l2 \ X ji2 , rjj 2 ) = rjj 2 h 0 2 (tj i2 ) exp{Xj l2 {3 2 }, t ji2 >0 
h 3 (tji 2 \tjii,Xji 3 ,Tjj 3 ') f/j3^03( tjl‘2 tjil) -^ 713 /^ 3 }) tji2 tjilj (17) 

Key features of this model, in relation to the proposed framework are: 

• Cluster-specific effects are represented via the (rjj i, rjj 2 , rjj 3 ), each of which is assigned 
an independent univariate parametric distribution (either a log-Normal or a Gamma). 
As such, the model does not permit the characterization of covariation between the 
cluster-specific random effects. In contrast, the proposed methods provides analysts 
with two choices for the joint distribution of the V)’s: a parametric MVN or a non- 
parametric DPM-MVN. 

• There is no patient-specific term, analogous to the 7 ^ in the proposed model. As 
such a potentially important source of within-subject correlation between 7\ and T 2 
is not accounted for. 

• Similar to the propose methods, however, is that the baseline hazard function for 
h 3 () can be specified non-parametrically (via a spline) or parametrically (using the 
Weibull distribution). 

• Although not evident from the model specification, estimation of the shared frailty 
model is based on three separate fits of the three models. In contrast, because the 
proposed model considers several components of covariation (i.e. covariation among 
the V/s and the patient-specific Tp’s) we perform estimation/inference using single 
likelihood. Indeed for the shared frailty model to accommodate these components of 
covariation, a completely new framework for estimation/inference would need to be 
developed. 

• Estimation of the (rjji, rjj 2l rjj 3 ) proceeds using empirical Bayes (after estimation of the 
remaining components via an integrated likelihood). Uncertainty for these estimates 
are only available when their distributions are taken to be Gamma distributions. 
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C.1.2 The joint frailty model 

Two variations of a joint frailty model have been implemented in the frailtypack package. 
The first was developed for the analysis of a recurrent non-terminal event and a terminal 
event and specifies a single hazard function for each. Specifically, the model is given by: 

hi{tkn\u>i) = ayr 0 (t ki] ) exp{ Xj[ /3]}, for recurrent non-terminal event 
h 2 (ti 2 \wi) = exp{X^/3 2 }, for the terminal event 

u>i ~ Gamma(l/0, 1/9). (18) 

where ay is a common subject-specihc frailty representing unobserved covariates that impact 
both events. We note that this is specification is similar to the model proposed by Liu et ah 
(2004). 

The second joint frailty model implemented in the frailtypack package is for model¬ 
ing two clustered survival outcomes. Specifically, the model posits that the event-specific 
hazard functions for the j th cluster are structured as follows: 

= ho^tjn) expiX^Pi + rjj}, for any event 
h(tji 2 \ r )j) = ho 2 (tji 2 ) exp{A% 2 /3 2 + m^}, for the terminal event 

r]j ~ Normal(0, a 2 ) (19) 

In relation to the context we consider, the central limitation of these model is that they 
only consider a single level of the two-level hierarchy inherent to cluster-correlated semi- 
competing risks data. Specifically, as applied and described in the fraitypack package, 
the first model only considers patient-specific effects while the second model only considers 
cluster-specific effects. As such neither model would be appropriate for our motivating 
application since (i) ignoring cluster-specific effects means that one cannot address several 
of our key scientific questions and (ii) ignoring patient-level effects can result in substantial 
bias (see the simulation studies in Section 5 of the main manuscript). 

We also note that a second limitation is that model (19) does not consider the transition 
from the non-terminal event to the terminal event; that is there is no analogue for () in 
the model. This represents a limitation in the sense that information readily available in 
the data is ignored. In the motivating application in the main manuscript, for example, the 
fact that the time of death following readmission within 90 days is known for 608 (11.5%) 
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patients is ignored. Finally, although model (18) does permit a patient to transition from 
the non-terminal state to the terminal state, this transition is assumed to occur at the 
same rate at which a patient who is in the initial state transitions directly into the terminal 
state; that is, in contrast to the proposed model that distinguishes h 2 () from /i 3 (), model 
(18) only has a single hazard for the terminal event. 

C.2 Gorfine and Hsu (2011) 

Gorfine and Hsu (2011) explicitly consider the related but distinct problem of analyzing 
cluster correlated competing risk data for which Tj and T 2 are both terminal events (i.e. 
death due to two causes). Towards analyzing such data, they propose the following hierar¬ 
chical model: 


hitjnlXj^ej^tjii)) = h 01 (tjii) exp{Xf i (3i +for cause 1 
h(t ji2 \X jU e j2 (t ji2 )) = h 02 (t ji2 ) exp{Xj ^2 + e j2 (t ji2 )}, for cause 2 (20) 

to describe the risk of transitioning into one of the two terminal states for the i th patient in 
the j th cluster. As part of their development, Gorfine and Hsu (2011) provide a framework 
within which the distribution of the cluster-specific Cj g (t ) terms can be flexibly specified. 

While this flexibility is very appealing, direct application of this model to our motivating 
application would be subject to a number of limitations mainly because the method was 
not designed for the cluster-correlated semi-competing risks setting. Specifically, 

• Similar to the joint frailty model given by (19), the application of model (20) means 
that one would ignore information in the data on the transition from the non-terminal 
event to the terminal event; that is, there is not analogue of /i 3 (). 

• Although model (20) includes cluster-specific random effects, it does not include 
specification of patient-specific terms analogous to 7 ^ in the proposed model. As is 
clear from the simulations presented in Section 5 of the main manuscript, ignoring this 
component of variation can lead to substantial bias in estimation and poor inferential 
properties in the cluster-correlated semi-competing risks setting. 
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D Simulation Results 


In order to supplement the results from simulation studies, we provide estimated percent 
bias, coverage probability, and average relative width of 95% credible/confidence inter¬ 
vals for /3 X , /3 2 , 0 3 , and 6 for our four proposed models and four types of SF models of 
Liquet et al. (2012) in Table D.1-D.6. We also provide estimated transition-specific base¬ 
line survival functions for the models under simulation scenarios 2,3, and 5 in Figure D.l. 
Note that since results from SF models are almost identical between models that adopt 
the independent gamma distributions for cluster-specific random effects and those that 
adopt the independent log-Normal distributions, we only present the results from SF mod¬ 
els with the gamma cluster-specific random effects in Figure D.l. We also present Table 
D.7 that augments Table 6 in the main manuscript by additionally presenting results for 
the Liquet et al. (2012)’s models that adopt independent log-Normal distributions for the 
cluster-specific random effects. 

The results presented in this section are generally consistent with the conclusions we 
drew in the main paper: contrary to the existing SF models, our proposed models yielded a 
small bias and coverage probability estimated to be close to the nominal 0.95 for regression 
parameters and 6 (except scenario 4 for which 6= 0); all four of the proposed models estimate 
the three baseline survival functions very well. 
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Tabic D.l: Estimated percent bias and coverage probability for ( 3 \ and 6 for six analyses described in Section 5.2, across six 
simulation scenarios given in Table 3. Throughout values are based on results from i?=500 simulated datasets. 


Percent Bias Coverage Probability 


Scenario 


True 

Weibull 

Weibull 

Weibull 

Weibull 

PEM 

PEM 

Spline 

Spline 

Weibull 

Weibull 

Weibull 

Weibull 

PEM 

PEM 

Spline 

Spline 



value 

-MVN 

-DPM 

-SF g a 

-SF ^ 

-MVN 

-DPM 

-SF g 

-SF CM 

-MVN 

-DPM 

-SF g 

-SF^ 

-MVN 

-DPM 

-SFg 

-SF^ 


fin 

0.50 

0.1 

0.2 

-19.8 

-21.3 

0.4 

0.4 

-21.0 

-20.8 

0.96 

0.96 

0.01 

0.01 

0.95 

0.96 

0.00 

0.00 

1 

fil2 

0.80 

0.2 

0.3 

-19.7 

-21.3 

0.5 

0.4 

-21.0 

-20.8 

0.95 

0.95 

0.00 

0.00 

0.96 

0.97 

0.00 

0.00 


fill 

-0.50 

0.3 

0.3 

-19.8 

-18.8 

0.3 

0.3 

-21.2 

-20.9 

0.97 

0.96 

0.31 

0.31 

0.96 

0.96 

0.25 

0.26 


e 

0.50 

1.0 

1.3 



1.4 

1.2 



0.95 

0.95 



0.93 

0.94 




fii i 

0.50 

-0.1 

-0.0 

-31.8 

-33.4 

0.1 

0.1 

-32.8 

-32.8 

0.94 

0.94 

0.00 

0.00 

0.94 

0.93 

0.00 

0.00 

2 

fil 2 

0.80 

0.1 

0.2 

-31.7 

-33.3 

0.4 

0.3 

-32.7 

-32.7 

0.97 

0.97 

0.00 

0.00 

0.94 

0.95 

0.00 

0.00 


fill 

-0.50 

1.2 

1.3 

-31.1 

-29.2 

1.1 

1.1 

-32.2 

-32.2 

0.94 

0.95 

0.05 

0.05 

0.94 

0.94 

0.04 

0.04 


e 

1.00 

0.4 

0.7 



0.7 

0.6 



0.94 

0.95 



0.94 

0.95 




fii i 

0.50 

0.3 

0.3 

-19.9 

-20.7 

0.7 

0.7 

-21.0 

-20.9 

0.94 

0.94 

0.00 

0.00 

0.93 

0.94 

0.00 

0.00 

3 

fin 

0.80 

0.4 

0.4 

-19.8 

-20.7 

0.8 

0.8 

-20.9 

-20.8 

0.94 

0.94 

0.00 

0.00 

0.94 

0.94 

0.00 

0.00 


fin 

-0.50 

0.4 

0.3 

-20.1 

-19.7 

0.5 

0.6 

-21.2 

-21.2 

0.96 

0.96 

0.31 

0.29 

0.95 

0.96 

0.27 

0.27 


e 

0.50 

2.0 

2.1 



3.2 

3.2 



0.96 

0.96 



0.93 

0.95 




fin 

0.50 

3.7 

3.7 

0.2 

-2.9 

4.7 

4.6 

0.3 

0.6 

0.87 

0.86 

0.96 

0.91 

0.81 

0.83 

0.96 

0.95 

4 

fin 

0.80 

3.6 

3.6 

-0.0 

-3.1 

4.5 

4.5 

0.1 

0.4 

0.80 

0.79 

0.95 

0.89 

0.69 

0.70 

0.95 

0.95 


fin 

-0.50 

4.0 

4.0 

0.2 

7.3 

4.8 

4.7 

0.2 

0.6 

0.93 

0.94 

0.94 

0.88 

0.93 

0.93 

0.93 

0.94 


e 

0.00 


















fin 

0.50 

-0.3 

0.1 

-20.3 

-24.6 

0.0 

0.3 

-21.1 

-20.9 

0.94 

0.95 

0.00 

0.00 

0.96 

0.96 

0.00 

0.00 

5 

fill 

0.80 

0.0 

0.3 

-20.0 

-24.6 

0.3 

0.6 

-20.9 

-20.7 

0.95 

0.95 

0.00 

0.00 

0.96 

0.96 

0.00 

0.00 


fill 

-0.50 

-0.2 

0.2 

-20.4 

-13.7 

-0.2 

0.2 

-21.3 

-21.1 

0.94 

0.94 

0.29 

0.26 

0.94 

0.94 

0.25 

0.26 


8 

0.50 

-0.2 

1.0 



0.4 

1.3 



0.95 

0.95 



0.95 

0.96 




fin 

0.50 

9.3 

9.4 

-22.1 

-23.6 

0.4 

0.3 

-25.9 

-25.1 

0.58 

0.57 

0.00 

0.00 

0.94 

0.94 

0.00 

0.00 

6 

fin 

0.80 

9.7 

9.8 

-22.0 

-23.5 

0.5 

0.5 

-25.8 

-25.0 

0.20 

0.20 

0.00 

0.00 

0.94 

0.95 

0.00 

0.00 


fin 

-0.50 

10.2 

10.2 

-21.6 

-18.2 

0.8 

0.7 

-26.1 

-24.9 

0.81 

0.80 

0.21 

0.21 

0.93 

0.94 

0.10 

0.10 


8 

0.50 

52.8 

53.0 



1.8 

1.7 



0.00 

0.00 



0.95 

0.96 




a The SF models that adopt the independent gamma distributions for cluster-specific random effects 
b The SF models that adopt the independent log-Normal distributions for cluster-specific random effects 












Table D.2: Estimated percent bias and coverage probability for ( 3 2 for six analyses described in Section 5.2, across six simulation 
scenarios given in Tabic 3. Throughout values are based on results from f?=500 simulated datasets. 

Percent Bias Coverage Probability 


Scenario 


True 

value 

Weibull 

-MVN 

Weibull 

-DPM 

Weibull 

-SFg“ 

Weibull 

-SF^ 6 

PEM 

-MVN 

PEM 

-DPM 

Spline 

-SFg 

Spline 

-sf £ ^ 

Weibull 

-MVN 

Weibull 

-DPM 

Weibull 

-SFg 

Weibull 

-SFg^ 

PEM 

-MVN 

PEM 

-DPM 

Spline 

-SFg 

Spline 

-SF CU 


021 

0.50 

-0.1 

-0.0 

-27.9 

-25.9 

-0.2 

-0.3 

-26.1 

-26.1 

0.93 

0.93 

0.00 

0.00 

0.94 

0.94 

0.00 

0.00 

1 

022 

0.80 

0.3 

0.4 

-27.9 

-25.6 

0.2 

0.1 

-25.8 

-25.7 

0.96 

0.96 

0.00 

0.00 

0.96 

0.96 

0.00 

0.00 


023 

-0.50 

1.3 

1.4 

-26.5 

-21.1 

0.9 

0.9 

-25.2 

-25.1 

0.94 

0.93 

0.34 

0.35 

0.95 

0.95 

0.34 

0.34 


021 

0.50 

-0.1 

0.1 

-39.2 

-39.2 

-0.2 

-0.2 

-39.9 

-39.8 

0.96 

0.96 

0.00 

0.00 

0.96 

0.96 

0.00 

0.00 

2 

022 

0.80 

0.2 

0.3 

-39.1 

-39.0 

0.0 

0.0 

-39.6 

-39.6 

0.96 

0.97 

0.00 

0.00 

0.96 

0.97 

0.00 

0.00 


023 

-0.50 

1.3 

1.4 

-38.5 

-38.3 

0.8 

0.8 

-39.2 

-39.2 

0.93 

0.94 

0.07 

0.07 

0.93 

0.93 

0.06 

0.06 


021 

0.50 

0.3 

0.3 

-27.4 

-24.9 

0.3 

0.3 

-25.7 

-25.6 

0.95 

0.95 

0.01 

0.01 

0.95 

0.96 

0.00 

0.00 

3 

022 

0.80 

0.5 

0.5 

-27.5 

-24.8 

0.5 

0.5 

-25.5 

-25.4 

0.95 

0.95 

0.00 

0.00 

0.94 

0.95 

0.00 

0.00 


023 

-0.50 

2.5 

2.5 

-24.9 

-23.2 

2.3 

2.5 

-24.2 

-24.1 

0.95 

0.95 

0.37 

0.38 

0.94 

0.94 

0.35 

0.36 


021 

0.50 

4.5 

4.5 

-4.9 

-4.4 

4.7 

4.7 

-0.6 

-0.6 

0.89 

0.89 

0.91 

0.90 

0.87 

0.88 

0.96 

0.96 

4 

022 

0.80 

4.6 

4.6 

-5.2 

-4.6 

4.8 

4.8 

-0.6 

-0.5 

0.78 

0.78 

0.88 

0.88 

0.76 

0.79 

0.92 

0.92 


023 

-0.50 

5.4 

5.4 

25.8 

18.1 

5.4 

5.5 

-0.1 

-0.0 

0.92 

0.92 

0.88 

0.88 

0.91 

0.93 

0.93 

0.92 


021 

0.50 

-0.3 

0.2 

-26.6 

-25.1 

-0.4 

-0.0 

-26.1 

-25.8 

0.94 

0.95 

0.00 

0.00 

0.94 

0.95 

0.00 

0.00 

5 

022 

0.80 

-0.4 

0.1 

-26.8 

-25.3 

-0.5 

-0.1 

-26.2 

-25.9 

0.96 

0.96 

0.00 

0.00 

0.96 

0.96 

0.00 

0.00 


023 

-0.50 

0.1 

0.5 

-27.2 

-24.9 

-0.4 

0.1 

-26.1 

-25.8 

0.95 

0.95 

0.33 

0.35 

0.95 

0.95 

0.30 

0.31 


021 

0.50 

9.2 

9.3 

-25.1 

-23.8 

0.1 

0.1 

-25.0 

-25.0 

0.70 

0.70 

0.01 

0.01 

0.95 

0.95 

0.00 

0.01 

6 

022 

0.80 

9.7 

9.7 

-25.4 

-23.8 

0.4 

0.3 

-24.8 

-24.8 

0.40 

0.39 

0.00 

0.00 

0.94 

0.95 

0.00 

0.00 


023 

-0.50 

10.4 

10.4 

-26.6 

-13.6 

0.8 

0.8 

-24.5 

-24.5 

0.85 

0.85 

0.47 

0.47 

0.94 

0.95 

0.36 

0.35 


a The SF models that adopt the independent gamma distributions for cluster-specific random effects 
b The SF models that adopt the independent log-Normal distributions for cluster-specific random effects 












Table D.3: Estimated percent bias and coverage probability for ( 3 3 for six analyses described in Section 5.2, across six simulation 
scenarios given in Tabic 3. Throughout values are based on results from i?=500 simulated datasets. 

Percent Bias Coverage Probability 


Scenario 


True 

value 

Weibull 

-MVN 

Weibull 

-DPM 

Weibull 

-SFg“ 

Weibull 

-SF^ 6 

PEM 

-MVN 

PEM 

-DPM 

Spline 

-SFg 

Spline 

-sf £ ^ 

Weibull 

-MVN 

Weibull 

-DPM 

Weibull 

-SFg 

Weibull 

-SFg^ 

PEM 

-MVN 

PEM 

-DPM 

Spline 

-SFg 

Spline 

-SF^ 


0i 31 

1.00 

0.4 

0.5 

-21.8 

-12.5 

0.7 

0.8 

-13.3 

-13.2 

0.94 

0.94 

0.08 

0.09 

0.94 

0.94 

0.06 

0.06 

1 

032 

1.00 

0.2 

0.3 

-20.5 

-9.0 

0.6 

0.6 

-9.7 

-9.7 

0.96 

0.96 

0.28 

0.32 

0.93 

0.94 

0.27 

0.28 


033 

-1.00 

0.2 

0.3 

44.3 

-12.5 

0.4 

0.4 

-13.4 

-13.3 

0.94 

0.94 

0.47 

0.53 

0.94 

0.94 

0.49 

0.50 


031 

1.00 

0.1 

0.3 

-24.1 

-25.3 

0.6 

0.6 

-24.8 

-24.7 

0.95 

0.95 

0.00 

0.00 

0.95 

0.95 

0.00 

0.00 

2 

032 

1.00 

0.4 

0.5 

-22.8 

-24.1 

0.8 

0.8 

-23.3 

-23.3 

0.94 

0.94 

0.00 

0.00 

0.95 

0.95 

0.00 

0.00 


033 

-1.00 

0.2 

0.4 

-23.8 

-22.6 

0.4 

0.4 

-24.6 

-24.5 

0.96 

0.96 

0.08 

0.07 

0.95 

0.95 

0.06 

0.07 


031 

1.00 

0.6 

0.6 

-21.4 

-14.4 

1.1 

1.1 

-12.6 

-12.5 

0.96 

0.96 

0.10 

0.14 

0.95 

0.95 

0.09 

0.09 

3 

032 

1.00 

0.7 

0.8 

-19.2 

-11.2 

1.2 

1.2 

-9.0 

-8.9 

0.95 

0.95 

0.35 

0.40 

0.94 

0.94 

0.37 

0.38 


033 

-1.00 

0.3 

0.3 

13.1 

-10.1 

0.5 

0.6 

-12.6 

-12.5 

0.95 

0.95 

0.52 

0.57 

0.93 

0.94 

0.54 

0.54 


031 

1.00 

3.4 

3.4 

-0.6 

11.9 

4.1 

3.9 

11.0 

11.3 

0.87 

0.88 

0.11 

0.12 

0.84 

0.86 

0.15 

0.15 

4 

032 

1.00 

3.4 

3.5 

5.6 

19.2 

4.2 

4.0 

18.3 

18.6 

0.88 

0.88 

0.00 

0.00 

0.83 

0.86 

0.00 

0.00 


033 

-1.00 

3.0 

3.0 

20.6 

11.6 

3.3 

3.5 

10.3 

10.7 

0.95 

0.95 

0.51 

0.57 

0.94 

0.94 

0.62 

0.61 


031 

1.00 

-0.0 

0.5 

-22.4 

-16.0 

0.4 

0.8 

-14.2 

-14.0 

0.97 

0.96 

0.05 

0.07 

0.96 

0.97 

0.04 

0.05 

5 

032 

1.00 

0.0 

0.6 

-20.2 

-12.4 

0.4 

0.9 

-10.4 

-10.1 

0.96 

0.95 

0.27 

0.31 

0.95 

0.95 

0.25 

0.26 


033 

-1.00 

0.0 

0.5 

36.3 

-11.4 

0.3 

0.6 

-14.1 

-13.9 

0.94 

0.94 

0.44 

0.49 

0.95 

0.94 

0.46 

0.47 


031 

1.00 

8.4 

8.5 

-28.5 

-28.5 

0.5 

0.5 

-30.4 

-30.3 

0.31 

0.30 

0.00 

0.00 

0.96 

0.95 

0.00 

0.00 

6 

032 

1.00 

8.9 

9.0 

-20.3 

-20.3 

0.6 

0.6 

-22.0 

-21.8 

0.28 

0.27 

0.00 

0.00 

0.95 

0.95 

0.00 

0.00 


033 

-1.00 

8.8 

8.9 

-27.9 

-27.8 

0.4 

0.4 

-30.2 

-30.1 

0.67 

0.67 

0.00 

0.00 

0.95 

0.95 

0.00 

0.00 


a The SF models that adopt the independent gamma distributions for cluster-specific random effects 
b The SF models that adopt the independent log-Normal distributions for cluster-specific random effects 
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Figure D.l: Estimated transition-specific baseline survival functions, So 5 (-)=exp(— H 0g (-)), 
for each six analyses described in Section 5 under simulation scenarios 2, 3 and 5. 
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Table D.4: Average relative width of 95% credible/confidence intervals for (3\ and 9, with 
the Weibull-MVN model taken as the referent, across six simulation scenarios given in 
Table 3. Throughout values are based on results from i?=500 simulated datasets. 


Scenario 


Weibull 

Weibull 

Weibull 

Weibull 

PEM 

PEM 

Spline 

Spline 



-MVN 

-DPM 

-SF g “ 

-SFV 

-MVN 

-DPM 

-SF e 

-SF £A f 


Pn 

1.00 

1.00 

0.81 

0.78 

1.02 

1.02 

0.81 

0.81 

1 

@12 

1.00 

1.00 

0.77 

0.74 

1.04 

1.04 

0.77 

0.77 


@13 

1.00 

1.00 

0.84 

0.81 

1.00 

1.01 

0.83 

0.84 


9 

1.00 

1.00 



1.10 

1.12 




@ii 

1.00 

1.00 

0.73 

0.70 

1.02 

1.02 

0.73 

0.73 

2 

P 12 

1.00 

1.00 

0.69 

0.66 

1.03 

1.04 

0.69 

0.69 


/ 5 l 3 

1.00 

1.00 

0.76 

0.72 

1.00 

1.00 

0.76 

0.76 


9 

1.00 

1.00 



1.12 

1.14 




P n 

1.00 

1.00 

0.81 

0.78 

1.02 

1.02 

0.81 

0.81 

3 

@12 

1.00 

1.00 

0.76 

0.74 

1.04 

1.04 

0.77 

0.77 


@13 

1.00 

1.00 

0.83 

0.80 

1.00 

1.01 

0.83 

0.83 


9 

1.00 

1.00 



1.10 

1.13 




@11 

1.00 

1.00 

0.95 

0.90 

1.02 

1.01 

0.96 

0.96 

4 

@12 

1.00 

1.00 

0.94 

0.88 

1.03 

1.03 

0.95 

0.95 


@13 

1.00 

1.00 

0.96 

0.90 

1.01 

1.01 

0.96 

0.96 


9 

1.00 

1.00 



1.09 

1.09 




@11 

1.00 

1.00 

0.81 

0.74 

1.02 

1.02 

0.81 

0.81 

5 

@12 

1.00 

1.00 

0.77 

0.70 

1.03 

1.03 

0.77 

0.77 


@13 

1.00 

1.00 

0.83 

0.76 

1.00 

1.00 

0.83 

0.83 


9 

1.00 

1.00 



1.09 

1.09 




@11 

1.00 

1.00 

0.74 

0.71 

0.94 

0.95 

0.73 

0.74 

6 

@12 

1.00 

1.00 

0.72 

0.69 

0.96 

0.97 

0.71 

0.72 


@13 

1.00 

1.00 

0.76 

0.72 

0.93 

0.93 

0.75 

0.75 


9 

1.00 

1.00 



0.89 

0.90 




“ The SF models that adopt the independent gamma distributions for cluster-specific random effects 
b The SF models that adopt the independent log-Normal distributions for cluster-specific random effects 










Table D.5: Average relative width of 95% credible/confidence intervals for fa, with the 
Weibull-MVN model taken as the referent, across six simulation scenarios given in Table 
3. Throughout values are based on results from i?=500 simulated datasets. 


Scenario 


Weibull 

-MVN 

Weibull 

-DPM 

Weibull 

-SF g a 

Weibull 

- SF ^ 

PEM 

-MVN 

PEM 

-DPM 

Spline 

- SF e 

Spline 

- SF £A , 


@21 

1.00 

1.00 

0.80 

0.82 

1.01 

1.01 

0.84 

0.84 

1 

@22 

1.00 

1.00 

0.75 

0.78 

1.02 

1.02 

0.79 

0.79 


@23 

1.00 

1.00 

0.83 

0.85 

1.00 

1.00 

0.87 

0.87 


@21 

1.00 

1.00 

0.77 

0.77 

1.01 

1.01 

0.78 

0.78 

2 

/?22 

1.00 

1.00 

0.72 

0.72 

1.02 

1.02 

0.73 

0.73 


@23 

1.00 

1.00 

0.80 

0.80 

1.00 

1.00 

0.80 

0.80 


@21 

1.00 

1.00 

0.80 

0.83 

1.01 

1.01 

0.83 

0.83 

3 

/?22 

1.00 

1.00 

0.75 

0.78 

1.02 

1.03 

0.79 

0.79 


@23 

1.00 

1.00 

0.82 

0.86 

1.00 

1.00 

0.86 

0.86 


@ 21 

1.00 

1.00 

0.90 

0.90 

1.01 

1.01 

0.96 

0.96 

4 

/?22 

1.00 

1.00 

0.88 

0.88 

1.01 

1.01 

0.95 

0.95 


@23 

1.00 

1.00 

0.91 

0.91 

1.00 

1.00 

0.97 

0.97 


@21 

1.00 

1.00 

0.81 

0.83 

1.01 

1.01 

0.84 

0.84 

5 

@22 

1.00 

1.00 

0.77 

0.78 

1.02 

1.02 

0.79 

0.79 


@23 

1.00 

1.00 

0.84 

0.86 

1.00 

1.00 

0.86 

0.86 


@21 

1.00 

1.00 

0.78 

0.79 

0.96 

0.96 

0.82 

0.82 

6 

@22 

1.00 

1.00 

0.75 

0.76 

0.97 

0.97 

0.80 

0.80 


@23 

1.00 

1.00 

0.79 

0.80 

0.95 

0.95 

0.84 

0.84 


“ The SF models that adopt the independent gamma distributions for cluster-specific random effects 
b The SF models that adopt the independent log-Normal distributions for cluster-specific random effects 
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Tabic D.6: Average relative width of 95% credible/confidence intervals for (3$, with the 
Weibull-MVN model taken as the referent, across six simulation scenarios given in Table 
3. Throughout values are based on results from i?=500 simulated datasets. 


Scenario 


Weibull 

-MVN 

Weibull 

-DPM 

Weibull 

-SF g a 

Weibull 

- SF m b 

PEM 

-MVN 

PEM 

-DPM 

Spline 

- SF e 

Spline 

- SF £A , 


@31 

1.00 

1.00 

0.72 

0.83 

1.01 

1.01 

0.83 

0.83 

1 

@32 

1.00 

1.00 

0.73 

0.83 

1.01 

1.02 

0.83 

0.83 


@33 

1.00 

1.00 

0.75 

0.86 

1.00 

1.00 

0.86 

0.86 


@31 

1.00 

1.00 

0.77 

0.75 

1.01 

1.01 

0.77 

0.77 

2 

/?32 

1.00 

1.00 

0.77 

0.75 

1.01 

1.02 

0.77 

0.77 


@33 

1.00 

1.00 

0.81 

0.79 

1.00 

1.00 

0.81 

0.81 


CO 

1.00 

1.00 

0.73 

0.81 

1.01 

1.01 

0.84 

0.84 

3 

@32 

1.00 

1.00 

0.74 

0.81 

1.01 

1.02 

0.84 

0.84 


@33 

1.00 

1.00 

0.76 

0.84 

1.00 

1.00 

0.87 

0.87 


co 

1.00 

1.00 

0.85 

0.97 

1.01 

1.01 

0.98 

0.98 

4 

@32 

1.00 

1.00 

0.87 

0.99 

1.01 

1.01 

1.00 

1.00 


@33 

1.00 

1.00 

0.85 

0.96 

1.00 

1.00 

0.97 

0.97 


@31 

1.00 

1.00 

0.73 

0.79 

1.01 

1.01 

0.83 

0.83 

5 

@32 

1.00 

1.00 

0.73 

0.80 

1.01 

1.01 

0.83 

0.84 


@33 

1.00 

1.00 

0.76 

0.82 

1.00 

1.00 

0.86 

0.86 


@31 

1.00 

1.00 

0.69 

0.69 

0.96 

0.96 

0.69 

0.69 

6 

/?32 

1.00 

1.00 

0.72 

0.72 

0.97 

0.97 

0.72 

0.72 


@33 

1.00 

1.00 

0.73 

0.73 

0.94 

0.94 

0.73 

0.73 


“ The SF models that adopt the independent gamma distributions for cluster-specific random effects 
b The SF models that adopt the independent log-Normal distributions for cluster-specific random effects 
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Table D.7: Mean squared error of prediction (xlCW 2 ) for cluster-specific random effects 
based on six analyses described in Section 5.2, across six data scenarios given in Table 3. 
Throughout values are based on results from i?=500 simulated datasets. 


Scenario 


Weibull 

-MVN 

Weibull 

-DPM 

Weibull 

-SFg“ 

%F 

Weibull 

-SF £A f fc 

%Ft 

PEM 

-MVN 

PEM 

-DPM 

Spline 

-SF s 

%F 

Spline 

-SF £A /- 

%F 


Vi 

5.25 

5.27 

6.40 


10355.30 


5.27 

5.27 

6.39 


10397.30 


1 

v 2 

7.66 

7.70 

8.70 

17.8 

11199.69 

6.6 

7.67 

7.72 

8.68 

0.2 

11131.71 

0.8 


v 3 

9.91 

9.95 

12.13 


11221.04 


9.91 

9.96 

12.11 


11214.38 



Vi 

6.36 

6.41 

8.10 


10535.62 


6.37 

6.41 

8.09 


10476.70 


2 

v 2 

8.76 

8.85 

10.23 

10.4 

11328.83 

7.8 

8.77 

8.86 

10.20 

0.0 

11208.54 

0.6 


V 3 

11.13 

11.19 

13.85 


11417.49 


11.13 

11.19 

13.91 


11449.00 



Vi 

5.03 

5.04 

6.27 


10398.49 


5.04 

5.04 

6.22 


10329.00 


3 

v 2 

6.34 

6.34 

8.28 

15.8 

11357.53 

9.2 

6.36 

6.36 

8.24 

0.0 

11348.16 

0.6 


V 3 

7.55 

7.49 

11.66 


10932.18 


7.57 

7.55 

11.69 


10895.77 



Vi 

3.84 

3.85 

4.99 


9798.48 


3.87 

3.87 

5.01 


9765.63 


4 

v 2 

6.25 

6.27 

7.19 

12.8 

11076.72 

12.4 

6.25 

6.27 

7.12 

0.4 

11102.66 

5.4 


V 3 

7.89 

7.90 

9.57 


10893.62 


7.90 

7.91 

9.52 


10886.70 



Vi 

6.95 

6.26 

10.87 


10005.24 


6.96 

6.27 

10.86 


9869.41 


5 

v 2 

11.52 

10.50 

14.95 

12.8 

11090.98 

78.4 

11.50 

10.52 

14.92 

0.2 

10976.78 

76.4 


V 3 

15.46 

14.66 

25.04 


11156.06 


15.46 

14.72 

24.94 


11073.46 



Vi 

5.05 

5.01 

6.34 


9670.25 


4.89 

4.85 

6.26 


9804.42 


6 

v 2 

7.58 

7.55 

8.60 

5.4 

11259.19 

9.2 

7.41 

7.39 

8.49 

1.4 

11272.07 

0.4 


V 3 

6.72 

6.65 

13.42 


10095.61 


6.44 

6.40 

13.70 


10233.22 



a The SF models that adopt the independent gamma distributions for cluster-specific random effects 
b The SF models that adopt the independent log-Normal distributions for cluster-specific random effects 
t % of times SF models yield at least one of Vj being — oo, resulting in MSEP being oo 
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E Application to Medicare data from New England 


In our main paper, posterior summaries for hazard ratio (HR) parameters for readmission, 
exp(/3i), from models for which a semi-Markov specification was adopted for h 03 (-) are 
presented. In Table E.1-E.4, we provide posterior summaries for HR parameters for death 
without readmission, exp(/3 2 ), and those for death following readmission, exp(/3 3 ), from 
both Markov and semi-Markov models. We also provide the posterior estimates of exp(/3i) 
from Markov models in Table E.5. 

From Table 4 (in the main paper) and Table E.1-E.5, we see the little difference in poste¬ 
rior estimates of HR parameters between Markov and semi-Markov models in this particular 
application. Therefore, our analyses in this document mainly focus on HR parameters for 
death (exp(/3 2 ) and exp(/3 3 )) under semi-Markov models. As seen in Table E.l and Table 
E.2, our proposed framework show how risk for death changes depending on whether or 
not a patient experiences the readmission. For example, whereas there is evidence of an 
increased risk of death for long hospital stay among individuals who have not been read¬ 
mitted (HR 1.10; 95% Cl 1.04, 1.18 in PEM-DPM), the same conclusion cannot be drawn 
for individuals who have been readmitted (HR 0.98; 95% Cl 0.87, 1.09 in PEM-DPM). 
In addition, the association between death and two of covariates (age and Charlson/Deyo 
score) is stronger in this magnitude (i.e. farther from zero) while the association between 
death and some other covariates (sex, source of entry to initial hospitalization, length of 
stay, discharged location, and whether or not patients underwent a procedure during the 
hospitalization) is weakened among patient who have been readmitted. In general, our 
analyses show evidence of increased risk for death for patients with male gender, older age, 
initially hospitalized via some route other than ER, higher comorbidity score, a procedure 
during the hospitalization, a discharge to a place other than home (without care). 
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Table E.l: Posterior medians (PM) and 95% credible intervals (Cl) for hazard ratio pa¬ 
rameters for death without readmission (exp(/3 2 )) from semi-competing risks data analyses 
based on semi-Markov models. 



Weibull-MVN 

PM (95%CI) 

Weibull-DPM 

PM (95%CI) 

PEM-MVN 

PM (95%CI) 

PEM-DPM 

PM (95%CI) 

Sex 

Male 

Female 

1.00 

0.69 (0.60, 0.79) 

1.00 

0.69 (0.61, 0.78) 

1.00 

0.75 (0.67, 0.83) 

1.00 

0.75 (0.67, 0.83) 

Age* 

1.09 (1.04, 1.13) 

1.09 (1.04, 1.14) 

1.07 (1.03, 1.11) 

1.07 (1.03, 1.11) 

Race 

White 

Non-white 

1.00 

0.93 (0.70, 1.22) 

1.00 

0.93 (0.70, 1.22) 

1.00 

0.94 (0.74, 1.17) 

1.00 

0.94 (0.75, 1.18) 

Source of entry to initial 

hospitalization 

Emergency room 

Other facility 

1.00 

1.61 (1.41, 1.85) 

1.00 

1.61 (1.41, 1.86) 

1.00 

1.50 (1.33, 1.70) 

1.00 

1.49 (1.32, 1.68) 

Charlson/Deyo score 

< 1 

> 1 

1.00 

1.40 (1.12, 1.71) 

1.00 

1.39 (1.13, 1.73) 

1.00 

1.26 (1.08, 1.50) 

1.00 

1.27 (1.06, 1.51) 

Procedure during hospitalization 

No 

Yes 

1.00 

0.09 (0.07, 0.12) 

1.00 

0.09 (0.07, 0.12) 

1.00 

0.13 (0.10, 0.16) 

1.00 

0.13 (0.10, 0.16) 

Length of stay** 

1.15 (1.07, 1.24) 

1.15 (1.06, 1.24) 

1.10 (1.04, 1.18) 

1.10 (1.04, 1.18) 

Care after discharge 

Home 

Home with care 

Hospice 

ICF/SNF 

Other 

1.00 

2.41 (2.00, 2.91) 
22.99 (18.08, 30.16) 
5.22 (4.29, 6.39) 
4.81 (3.48, 6.70) 

1.00 

2.45 (2.02, 2.94) 
23.71 (18.28, 31.20) 
5.33 (4.32, 6.45) 
4.93 (3.58, 6.84) 

1.00 

2.22 (1.85, 2.63) 
13.94 (11.22, 17.43) 
4.25 (3.57, 5.06) 
3.79 (2.91, 4.98) 

1.00 

2.21 (1.90, 2.61) 
13.85 (11.33, 17.08) 
4.25 (3.66, 5.01) 
3.81 (2.94, 4.91) 


* standardized so that 0 corresponds to an age of 77 years and so that one unit increment corresponds to 10 years 
** standardized so that 0 corresponds to 10 days and so that one unit increment corresponds to 7 days 
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Table E.2: Posterior medians (PM) and 95% credible intervals (Cl) for hazard ratio param¬ 
eters for death following readmission (exp(/3 3 )) from semi-competing risks data analyses 
based on semi-Markov models. 



Weibull-MVN 

PM (95%CI) 

Weibull-DPM 

PM (95%CI) 

PEM-MVN 

PM (95%CI) 

PEM-DPM 

PM (95%CI) 

Sex 

Male 

Female 

1.00 

0.81 (0.66, 1.00) 

1.00 

0.81 (0.66, 0.98) 

1.00 

0.84 (0.70, 1.00) 

1.00 

0.84 (0.71, 1.01) 

Age* 

1.10 (1.02, 1.19) 

1.10 (1.02, 1.19) 

1.09 (1.02, 1.17) 

1.09 (1.02, 1.17) 

Race 

White 

Non-white 

1.00 

1.15 (0.77, 1.67) 

1.00 

1.14 (0.79, 1.65) 

1.00 

1.12 (0.79, 1.54) 

1.00 

1.11 (0.78, 1.54) 

Source of entry to initial 

hospitalization 

Emergency room 

Other facility 

1.00 

1.54 (1.25, 1.90) 

1.00 

1.55 (1.25, 1.91) 

1.00 

1.42 (1.17, 1.72) 

1.00 

1.42 (1.16, 1.72) 

Charlson/Deyo score 

< 1 

> 1 

1.00 

1.51 (1.11, 2.06) 

1.00 

1.52 (1.11, 2.07) 

1.00 

1.41 (1.06, 1.85) 

1.00 

1.40 (1.05, 1.84) 

Procedure during hospitalization 

No 

Yes 

1.00 

0.21 (0.15, 0.29) 

1.00 

0.21 (0.15, 0.29) 

1.00 

0.28 (0.20, 0.39) 

1.00 

0.28 (0.21, 0.39) 

Length of stay** 

1.01 (0.89, 1.13) 

1.01 (0.89, 1.13) 

0.98 (0.88, 1.09) 

0.98 (0.87, 1.09) 

Care after discharge 

Home 

Home with care 

Hospice 

ICF/SNF 

Other 

1.00 

1.44 (1.13, 1.81) 
10.23 (4.66, 22.01) 
2.54 (1.87, 3.45) 
2.78 (1.64, 4.49) 

1.00 

1.44 (1.13, 1.82) 
10.43 (4.83, 22.33) 
2.57 (1.87, 3.46) 
2.72 (1.61, 4.44) 

1.00 

1.35 (1.08, 1.68) 
6.46 (3.33, 12.58) 
2.08 (1.52, 2.77) 
2.24 (1.40, 3.48) 

1.00 

1.34 (1.08, 1.65) 
6.35 (3.30, 12.29) 
2.07 (1.56, 2.76) 
2.25 (1.41, 3.43) 


* standardized so that 0 corresponds to an age of 77 years and so that one unit increment corresponds to 10 years 
** standardized so that 0 corresponds to 10 days and so that one unit increment corresponds to 7 days 
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Table E.3: Posterior medians (PM) and 95% credible intervals (Cl) for hazard ratio pa¬ 
rameters for death without readmission (exp(/3 2 )) from semi-competing risks data analyses 
based on Markov models. 



Weibull-MVN 

PM (95%CI) 

Weibull-DPM 

PM (95%CI) 

PEM-MVN 

PM (95%CI) 

PEM-DPM 

PM (95%CI) 

Sex 

Male 

Female 

1.00 

0.69 (0.60, 0.79) 

1.00 

0.69 (0.60, 0.79) 

1.00 

0.75 (0.67, 0.84) 

1.00 

0.75 (0.67, 0.83) 

Age* 

1.08 (1.04, 1.13) 

1.08 (1.04, 1.13) 

1.07 (1.03, 1.11) 

1.07 (1.03, 1.11) 

Race 

White 

Non-white 

1.00 

0.92 (0.69, 1.21) 

1.00 

0.92 (0.70, 1.22) 

1.00 

0.94 (0.75, 1.18) 

1.00 

0.94 (0.75, 1.16) 

Source of entry to initial 

hospitalization 

Emergency room 

Other facility 

1.00 

1.61 (1.42, 1.84) 

1.00 

1.62 (1.42, 1.86) 

1.00 

1.49 (1.32, 1.67) 

1.00 

1.49 (1.33, 1.69) 

Charlson/Deyo score 

< 1 

> 1 

1.00 

1.40 (1.12, 1.73) 

1.00 

1.39 (1.12, 1.72) 

1.00 

1.26 (1.06, 1.49) 

1.00 

1.27 (1.06, 1.51) 

Procedure during hospitalization 

No 

Yes 

1.00 

0.09 (0.07, 0.12) 

1.00 

0.09 (0.07, 0.12) 

1.00 

0.13 (0.10, 0.16) 

1.00 

0.13 (0.11, 0.17) 

Length of stay** 

1.15 (1.07, 1.23) 

1.15 (1.07, 1.23) 

1.10 (1.04, 1.17) 

1.10 (1.04, 1.17) 

Care after discharge 

Home 

Home with care 

Hospice 

ICF/SNF 

Other 

1.00 

2.44 (2.04, 2.91) 
23.52 (17.98, 30.13) 
5.30 (4.36, 6.43) 
4.88 (3.59, 6.72) 

1.00 

2.42 (2.02, 2.93) 
23.55 (18.05, 30.53) 
5.29 (4.38, 6.47) 
4.87 (3.54, 6.73) 

1.00 

2.20 (1.86, 2.62) 
13.72 (11.22, 17.49) 
4.23 (3.61, 5.16) 
3.78 (2.93, 4.99) 

1.00 

2.20 (1.90, 2.61) 
13.78 (11.10, 17.1) 
4.25 (3.62, 5.02) 
3.82 (2.92, 4.97) 


* standardized so that 0 corresponds to an age of 77 years and so that one unit increment corresponds to 10 years 
** standardized so that 0 corresponds to 10 days and so that one unit increment corresponds to 7 days 


35 





Table E.4: Posterior medians (PM) and 95% credible intervals (Cl) for hazard ratio param¬ 
eters for death following readmission (exp(/3 3 )) from semi-competing risks data analyses 
based on Markov models. 



Weibull-MVN 

PM (95%CI) 

Weibull-DPM 

PM (95%CI) 

PEM-MVN 

PM (95%CI) 

PEM-DPM 

PM (95%CI) 

Sex 

Male 

Female 

1.00 

0.81 (0.66, 0.98) 

1.00 

0.81 (0.67, 0.99) 

1.00 

0.84 (0.71, 1.01) 

1.00 

0.85 (0.71, 1.02) 

Age* 

1.11 (1.02, 1.19) 

1.10 (1.03, 1.20) 

1.10 (1.03, 1.18) 

1.09 (1.02, 1.18) 

Race 

White 

Non-white 

1.00 

1.14 (0.78, 1.64) 

1.00 

1.15 (0.79, 1.67) 

1.00 

1.13 (0.81, 1.55) 

1.00 

1.14 (0.79, 1.58) 

Source of entry to initial 

hospitalization 

Emergency room 

Other facility 

1.00 

1.58 (1.28, 1.97) 

1.00 

1.58 (1.28, 1.97) 

1.00 

1.44 (1.18, 1.75) 

1.00 

1.46 (1.21, 1.77) 

Charlson/Deyo score 

< 1 

> 1 

1.00 

1.53 (1.12, 2.11) 

1.00 

1.53 (1.12, 2.11) 

1.00 

1.40 (1.02, 1.84) 

1.00 

1.40 (1.06, 1.86) 

Procedure during hospitalization 

No 

Yes 

1.00 

0.20 (0.14, 0.28) 

1.00 

0.20 (0.14, 0.28) 

1.00 

0.27 (0.19, 0.37) 

1.00 

0.27 (0.19, 0.36) 

Length of stay** 

1.00 (0.89, 1.13) 

1.01 (0.89, 1.13) 

0.98 (0.88, 1.09) 

0.98 (0.88, 1.08) 

Care after discharge 

Home 

Home with care 

Hospice 

ICF/SNF 

Other 

1.00 

1.44 (1.15, 1.82) 
11.81 (5.18, 25.66) 
2.70 (1.96, 3.68) 
2.92 (1.74, 4.77) 

1.00 

1.44 (1.13, 1.81) 
11.6 (5.08, 24.49) 
2.69 (1.99, 3.61) 
2.89 (1.74, 4.68) 

1.00 

1.32 (1.06, 1.63) 
6.95 (3.49, 12.75) 

2.12 (1.59, 2.81) 

2.32 (1.46, 3.65) 

1.00 

1.33 (1.07, 1.66) 
6.79 (3.24, 13.28) 
2.17 (1.63, 2.87) 

2.36 (1.47, 3.67) 


* standardized so that 0 corresponds to an age of 77 years and so that one unit increment corresponds to 10 years 
** standardized so that 0 corresponds to 10 days and so that one unit increment corresponds to 7 days 
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Table E.5: Posterior medians (PM) and 95% credible intervals (Cl) for hazard ratio param¬ 
eters for readmission (exp(/3i)) from semi-competing risks data analyses based on Markov 
models. 



Weibull-MVN 

PM (95%CI) 

Weibull-DPM 

PM (95%CI) 

PEM-MVN 

PM (95%CI) 

PEM-DPM 

PM (95%CI) 

Sex 

Male 

Female 

1.00 

0.79 (0.70, 0.91) 

1.00 

0.80 (0.70, 0.91) 

1.00 

0.85 (0.76, 0.95) 

1.00 

0.85 (0.76, 0.95) 

Age* 

0.90 (0.86, 0.95) 

0.90 (0.86, 0.95) 

0.91 (0.87, 0.94) 

0.91 (0.87, 0.95) 

Race 

White 

Non-white 

1.00 

1.11 (0.86, 1.45) 

1.00 

1.11 (0.86, 1.44) 

1.00 

1.12 (0.89, 1.40) 

1.00 

1.11 (0.89, 1.38) 

Source of entry to initial 

hospitalization 

Emergency room 

Other facility 

1.00 

1.18 (1.03, 1.35) 

1.00 

1.19 (1.03, 1.36) 

1.00 

1.12 (0.99, 1.26) 

1.00 

1.12 (0.99, 1.27) 

Charlson/Deyo score 

< 1 

> 1 

1.00 

1.49 (1.19, 1.84) 

1.00 

1.50 (1.19, 1.85) 

1.00 

1.40 (1.15, 1.68) 

1.00 

1.39 (1.15, 1.68) 

Procedure during hospitalization 

No 

Yes 

1.00 

0.45 (0.37, 0.53) 

1.00 

0.45 (0.37, 0.53) 

1.00 

0.57 (0.49, 0.66) 

1.00 

0.57 (0.48, 0.66) 

Length of stay** 

1.15 (1.07, 1.23) 

1.15 (1.07, 1.23) 

1.12 (1.05, 1.19) 

1.12 (1.05, 1.19) 

Care after discharge 

Home 

Home with care 

Hospice 

ICF/SNF 

Other 

1.00 

0.95 (0.82, 1.11) 
0.39 (0.23, 0.62) 
0.88 (0.73, 1.06) 
1.05 (0.77, 1.43) 

1.00 

0.95 (0.82, 1.11) 
0.38 (0.22, 0.64) 
0.88 (0.73, 1.07) 
1.04 (0.77, 1.45) 

1.00 

0.89 (0.78, 1.02) 

0.27 (0.16, 0.42) 

0.76 (0.63, 0.90) 

0.89 (0.68, 1.18) 

1.00 

0.89 (0.78, 1.01) 
0.27 (0.16, 0.43) 
0.76 (0.64, 0.90) 
0.89 (0.67, 1.18) 


* standardized so that 0 corresponds to an age of 77 years and so that one unit increment corresponds to 10 years 
** standardized so that 0 corresponds to 10 days and so that one unit increment corresponds to 7 days 
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For our proposed models, we assess the convergence of our MCMC scheme by evaluating 
the potential scale reduction factor (PSRF) of Gelman et al. (2013). The potential problem 
with PSRF is that it has not converged but happens to be close to 1 by chance even though 
the PSRF is actually fluctuating. Therefore, for each parameter, the PSRF was calculated 
at several points in time with the first half discarded as burn-in. Then, we summarize 
the results using mean, maximum, and minimum value of PSRF for all model parameters 
at different iterations. The results are shown in Figure E.l. As the number of MCMC 
iterations increases, the mean PSRF converges toward 1 and the maximum of PSRF is less 
than 1.05 indicating that all model parameters have converged well. 
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Figure E.l: The mean, maximum, and minimum value of the potential scale reduction 
factor (PSRF) of all model parameters from the analysis of Medicare data. 
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